Structured penalties for functional linear models---partially empirical eigenvectors for regression

One of the challenges with functional data is incorporating spatial structure, or local correlation, into the analysis. This structure is inherent in the output from an increasing number of biomedical technologies, and a functional linear model is often used to estimate the relationship between the predictor functions and scalar responses. Common approaches to the ill-posed problem of estimating a coefficient function typically involve two stages: regularization and estimation. Regularization is usually done via dimension reduction, projecting onto a predefined span of basis functions or a reduced set of eigenvectors (principal components). In contrast, we present a unified approach that directly incorporates spatial structure into the estimation process by exploiting the joint eigenproperties of the predictors and a linear penalty operator. In this sense, the components in the regression are `partially empirical' and the framework is provided by the generalized singular value decomposition (GSVD). The GSVD clarifies the penalized estimation process and informs the choice of penalty by making explicit the joint influence of the penalty and predictors on the bias, variance, and performance of the estimated coefficient function. Laboratory spectroscopy data and simulations are used to illustrate the concepts.


Introduction
The coefficient function, β, in a functional linear model (fLM) represents the linear relationship between responses, y, and predictors, x, either of which may appear as a function. We consider the special case of scalar-on-function regression, formally written as y = I x(t)β(t) dt + ǫ, where x is a random function, square integrable on a closed interval I ⊂ Ê, and ǫ a vector of random i.i.d. meanzero errors. In many instances, one has an approximate idea about the informative structure of the predictors, such as the extent to which they are smooth, oscillatory, peaked, etc. Here we focus on analytical framework for incorporating such information into the estimation of β.
The analysis of data in this context involves a set of n responses {y i } n i=1 corresponding to a set of predictor curves {x i } n i=1 , each arising as a discretized sampling of an idealized function; i.e., x i ≡ (x i (t 1 ), ..., x i (t p )), for some, t 1 < · · · < t p , of I. In particular, the concept of geometric or spatial structure implies an order relation among the index parameter values. We assume the predictor functions have been sampled equally and densely enough to capture geometric structure of the type typically attributed to functions in (subspaces of) L 2 (I). For this, it will be assumed that p > n although this condition is not necessary for our discussion.
Several methods for estimating β are based on the eigenfunctions associated with the autocovariance operator defined by the predictors [16,32]. These eigenfunctions provide an empirical basis for representing the estimate and are the basis for the usual ordinary least-squares and principal-component estimates in multivariate analysis. The book by Ramsay and Silverman [38] summarize a variety of estimation methods that involve some combination of the empirical eigenfunctions and smoothing, using B-splines or other technique, but none of these methods provide an analytically tractable way to incorporate presumed structure directly into the estimation process. The approach presented here achieves this by way of a penalty operator, L, defined on the space of predictor functions.
The joint influence of the penalty and predictors on the estimated coefficient function is made explicit by way of the generalized singular value decomposition (GSVD) for a matrix pair. Just as the ordinary SVD provides the ingredients for an ordinary least squares estimate (in terms of the empirical basis), the GSVD provides a natural way to express a penalized least-squares estimate in terms of a basis derived from both the penalty and the predictors. We describe this in terms of the n × p matrix of sampled predictors, X, and an m × p discretized penalty operator, L. The general formulation is familiar as we consider estimates of β that arise from a squared-error loss with quadratic penalty:β α,L = arg min β {||y − Xβ|| 2 Ê n + α||Lβ|| 2 What distinguishes our presentation from others using this formulation is an emphasis on the joint spectral properties of the pair (X, L), as arise from the GSVD. We investigate the analytical role played by L in imposing structure on the estimate and focus on how the structure of L's leastdominant singular vectors should be commensurate with the informative structure of β.
In a Bayesian view, one may think of L as implementing a prior that favors a coefficient function lying near a particular subspace; this subspace is determined jointly by X and L. We note, however, that informative priors must come from somewhere and while they may come from expectations regarding smoothness, other information often exists-including pilot data, scientific knowledge or laboratory and instrumental properties. Our presentation aims to elucidate the role of L in providing a flexible means of implementing informative priors, regardless of their origin.
The general concept of incorporating "structural information" into regularized estimation for functional and image data is well established [2,12,36]. Methods for penalized regression have adopted this by constraining high-dimensional problems in various "structured" ways (sometimes with use of an L 1 norm): locally-constant structure [49,46], spatial smoothness [20], correlationbased constraints [52], and network-dependence structure described via a graph [26]. These general penalties have been motivated by a variety of heuristics: Huang et al. [24] refer to the seconddifference penalty as an "intuitive choice"; Hastie et al. [20] refer to a "structured penalty matrix [which] imposes smoothness with regard to an underlying space, time or frequency domain"; Tibshirani and Taylor [50] note that the rows of L should "reflect some believed structure or geometry in the signal"; and the penalties of Slawski et al. [46] aim to capture "a priori association structure of the features in more generality than the fused lasso." The most common penalty is a (discretized) derivative operator, motivated by the heuristic of penalizing roughness (see [21,38]). Our perspective on this is more analytical: since the eigenfunctions of the second-derivative operator L = D 2 (with zero boundary conditions on [0, 1]) are of the form ϕ(t) = sin(kπt), with eigenvalues k 2 π 2 (k = 1, 2, ..), L implements the assumption that the coefficient function is well represented by low-frequency trigonometric functions. This is in contrast to ridge regression (L = I) which imposes no geometric structure. Although not typically viewed this way, the choice of L = D 2 , or any differential operator, implies a favored basis for expansion of the estimate.
A purely empirical basis comprised of a few dominant right singular vectors of X is a common and theoretically tractable choice. This is the essence of principal component regression (PCR) and these vectors also form the basis for a ridge estimate. Although this empirical basis does not technically impose local spatial structure (no order relation among the index parameter values is used), it may be justified by arguing that a few principal component vectors capture the "greatest part" of a set of predictors [17]. Properties of this approach for signal regression is the focus of [7] and [16]. The functional data analysis framework of Ramsay and Silverman [38] provides two formulations of PCR. One in which the predictor curves are themselves smoothed prior to construction of principal components (chap. 8) and another that incorporates a roughness penalty into the construction of principal components (chap. 9), as originally proposed in [45]. In a related presentation on signal regression, Marx and Eilers [30] proposed a penalized B-spline approach in which predictors are transformed using a basis external to the problem (B-splines) and the estimated coefficient function is derived using the transform coefficients. Combining ideas from [30] and [21], the "smooth principal components" method of [8] projects predictors onto the dominant eigenfunctions to obtain an estimate then uses B-splines in a procedure that smooths the estimate. Reiss and Ogden [40] provide a thorough study on several of these methods and propose modifications that include two versions of PCR using B-splines and second-derivative penalties: FPCR C applies the penalty to the construction of the principal components (cf. [45]), while FPCR R incorporates the penalty into the regression (cf. [38]).
In the context of nonparametric regression (X = I) the formulation (1) plays a dominant role for smoothing [54]. Related to this, Heckman and Ramsay [22] proposed a differential equations modelbased estimate of a function µ whose properties are determined by a linear differential operator chosen from a parameterized family of differential equations, Lµ = 0. In this context, however, the GSVD is irrelevant since X does not appear and the role of L is relatively transparent.
Algebraic details on the GSVD as it relates to penalized least-squares are given in section 3 with analytic expressions for various properties of the estimation process are described in section 3.2. Intuitively, smaller bias is obtained by an informed choice of L (the goal being small Lβ). The affect of such a choice on the variance is described analytically. Section 4 describes several classes of structured penalties including two previously-proposed special cases that were justified by numerical simulations. The targeted penalties of subsection 4.2 are studied in more detail in section 5 including an analysis of the mean squared error for a family of penalized estimates which encompasses the ridge, principal-component and James-Stein estimates.
The assumptions on L here are increasingly restrictive to the point where the estimates are only minor extensions of these well-studied estimates. The goal, however, is to analytically describe the substantial gains achievable by even mild extensions of these established methods.
In applications the selection of the tuning parameter, α in (1), is important and so Section 6 describes our application of REML-based estimation for this. Numerical illustrations are provided in section 7: the simulation in subsection 7.1 is motivated by Reiss and Ogden's study of fLMs [40]; 7.2 presents a simulation using experimentally-derived Raman spectroscopy curves in which the "true" β has naturally-occurring (laboratory) structure; and section 7.3 presents an application based on experimentally collected spectroscopy curves representing varied biochemical (nanoparticle) concentrations. An appendix looks at the simulation studied by Hall and Horowitz [16]. We begin in section 2 with a brief setup for notation and an introductory example. Note that for any L = I, the estimated β is not given in terms of the ordinary empirical singular vectors (of X), but rather in terms of a "partially empirical" basis arising from a simultaneous diagonalization of X ′ X and L ′ L via the GSVD. Hence, for brevity, we refer toβ α,L as a PEER (partially empirical eigenvector for regression) estimate whenever L = I.

Background and simple example
Let β represent a linear functional on L 2 (I) defining a linear relationship y = I x(t)β(t) dt + ǫ (observed with error, ǫ) between a response, y, and random predictor function, x ∈ L 2 (I). We assume a set of n scalar responses {y i } n i=1 corresponding to the set of n predictors, {x i } n i=1 , each discretely sampled at common locations in I. Denote by X the n × p matrix whose ith row is a p-dimensional vector, x i , of discretely sampled functions, and columns that are centered to have mean 0. The notation ·, · will be used to denote the inner product on either L 2 (I) or Ê p , depending on the context. The empirical covariance operator is K = 1 n X ′ X, but for functional predictors, typically p > n or else K is ill-conditioned or rank deficient. In this case, there are either infinitely many least-squares solutions,β ≡ arg min β ||y − Xβ|| 2 , or else any such solution is highly unstable and of little use. The least-squares solution having minimum norm is unique, however, and it can be obtained directly by the singular value decomposition (SVD): X = U DV ′ where the left and right singular vectors, u k and v k , are the columns of U and V , respectively, and D = [D 1 0], where D 1 = diag{σ k } n k=1 , typically ordered as σ 1 ≥ σ 2 ≥ ... ≥ σ r > 0 (r = rank(X)). In terms of the SVD of X, the minimum- The orthogonal vectors that form the columns of V are the eigenvectors of X ′ X and are sometimes referred to as a Karhunen-Loève (K-L) basis for the row space of X.
The solutionβ + is Marquardt's generalized inverse estimator whose properties are discussed in [29]. For functional data,β + is an unstable, meaningless solution. One obvious fix is to truncate the sum to d < r terms so that {σ k } d k=1 is bounded away from zero. This leads to the truncated singular value or principal component regression (PCR) estimate: and subsequently, we use the notation A d ≡ col[a 1 , ..., a d ] to denote the first d columns of a matrix A.
When L = I, the minimizer in (1) is the ridge penalty due to A. E. Hoerl [23] or,β α, The factor F α acts to counterweight, rather than truncate, the terms 1 σ k as they get large. This is one of many possible filter factors which address problems of ill-determined rank (for more, see [12,19,33]). Weighted (or generalized) ridge regression replaces L = I with a diagonal matrix whose entries downweight those terms corresponding to the most variation [23]. Other "generalized ridge" estimates replace L = I by a discretized second-derivative operator, L = D 2 . Indeed, the Tikhonov-Phillips form of regularization (1) has a long history in the context of differential equations [51,36] and image analysis [15,33] with emphasis on numerical stability. In a linear model context, the smoothing imposed by this penalty was mentioned by Hastie and Mallows [21], discussed in Ramsay and Silverman [38] and used (on a the space of spline-transform coefficients) by Marx and Eilers [31], among others. The following simple example illustrates basic behavior for some of these penalties alongside an idealized PEER penalty.

A simple example
We consider a set of n = 50 bumpy predictor curves {x i } discretely sampled at p = 250 locations, as displayed in gray in the last panel of Figure 1. The true coefficient function, β, is displayed in black in this same panel. The responses are defined as y i = x i , β + ǫ i (ǫ i normal, uncorrelated mean-zero errors), and hence depend on the amplitudes of β's three bumps centered at locations t = 45, 110, 210.
A detailed simulation with complete results are provided in section 7.1. Here we simply illustrate the estimation process for L = I, as in (2), in comparison with L = D 2 and an idealized PEER penalty. The latter is constructed using a visual inspection of the predictors and lightly penalizes the subspace spanned by such structure, specifically, bumps centered at all visible locations (approximately t = 15, 45, 80, 110, 160, 210, 240).
The first five panels serve to emphasize the role played by the structure of basis vectors that comprise the series expansion in (2) (in terms of ordinary singular vectors) versus the analogous expansion (see (7)) in terms of generalized singular vectors. In particular, Figure 1 shows several partial sums of (7) for these three penalties. The ridge process (gray) is, naturally, dominated by the right singular vectors of X which become increasingly noisy in successive partial sums. The second-derivative penalized estimate (dashed) is dominated by low-frequency structure, while the targeted PEER estimate converges quickly to the informative features.
In this toy example, visual structure (spatial location) is used to define a regularization process that easily outperforms uninformed methods of penalization. Less visual examples where the penalty is defined by a set of laboratory-derived structure (in Raman spectroscopy curves) is given in sections 7.2 and 7.3; see Figure 2. In that setting, and in general, the role played by L is appropriately viewed in terms of a preferred subspace in Ê p determined by its singular vectors. Algebraic details about how structure in the estimation process is determined jointly by X and L = I are described next.

Penalized least squares and the GSVD
Of the many methods for estimating a coefficient function discussed in the Introduction, nearly are all aimed at imposing geometric or "functional" structure into the process via the use of basis functions in some manner. An alternative to choosing a basis outright is to exploit the structure imposed by an informed choice of penalty operator. The basis, determined by a pair (X, L), can be tailored toward structure of interest by the choice of L. When this is carried out in the leastsquares setting of (1), the algebraic properties of the GSVD explicitly reveal how the structure of the estimate is inherited from the spectral properties of (X, L).

The GSVD
For a given linear penalty L and parameter α > 0, the estimate in (1) takes the form This cannot be expressed using the singular vectors of X alone, but the generalized singular value decomposition of the pair (X, L) provides a tractable and interpretable series expansion. The GSVD appears in the literature in a variety of forms and notational conventions. Here we provide the necessary notation and properties of the GSVD for our purposes (see, e.g., [19]) but refer to [4,13,35] for a complete discussion and details about its computation. See also the comments of Bingham and Larntz [3]. Assume X is an n×p matrix (n ≤ p) of rank n, L is an m×p matrix (m ≤ p) of rank m. We also assume that n ≤ m ≤ p ≤ m + n, and the rank of the (n + m)× p matrix Z := [X ′ L ′ ] ′ is p. A unique solution is guaranteed if the null spaces of X and L intersect trivially: Null(L) ∩ Null(X) = {0}. This is not necessary for implementation, but it is natural in our applications and simplifies the notation. In addition, the condition p > n is not required, but rather than present notation for multiple cases, this will be assumed.
Given X and L, the following matrices exist and form the decomposition below: an n × n matrix U and an m × m matrix V , each with orthonormal columns, U ′ U = I, V ′ V = I; diagonal matrices S (n × n) and M (m × m); and a nonsingular p × p matrix W such that Here, S and M are of the form S = These matrices satisfy Denote the columns of U , V and W by u k , v k and w k , respectively. In spite of the notation, the generalized singular vectors u k and v k are not the same as the ordinary singular vectors of X in Section 2. They are the same when L = I, although their ordering is reversed; in that case, the ordinary singular values correspond to γ k := σ k /µ k for µ k > 0. By the convention used for ordering the GS values and vectors, the last few columns of W span the subspace Null(L) (or, if Null(L) is empty, they correspond to the smallest GS values, µ k ). We set d = dim(Null(L)) and note that µ k = 0 for k > n − d.
Now, equation (6) and some algebra gives ( A consequence of the ordering adopted for the GS values and vectors is that the first p − n columns of W play no role in this solution; see equation (4). So we can replace W by the p × n matrix, W n , consisting of the last n columns of W (corresponding to the indices p − n + 1 to p). We index the columns of W n as w 1 ,...,w n which is consistent with the indexing established in (5) for the singular values in S and M . Therefore, the L-penalized estimate can be expressed as a series in terms of GS values as This GSV expansion corresponds to a new basis for the estimation process: the estimate is expressed in terms of GS vectors {w k } determined jointly by X and L; cf. the ridge estimate in (2). For more compact notation used later, define the (n−d) For brevity, set o := n − d and let A o denote the first o columns of a matrix A. Also, denote by A ø the last d columns of A. In particular, the range of W ø is Null(L). Using this notation, (7) may be written concisely asβ where . In summary, the utility of a penalty L depends on whether the true coefficient function shares structural properties with this GSVD basis, {w k } n k=1 . With regard to this, the importance of the parameter α may be reduced by a judicious choice of L since the terms in (7) corresponding to the vectors {w k : µ k = 0} are independent of the parameter α [53].
As we'll see, bias enters the estimate to the extent that the vectors {w k : µ k = 0} appear in the expansion (7). The portion ofβ α,L that extends beyond the subspace Null(L) is constrained by a sphere (of radius determined by α); this portion corresponds to bias. Hence, L may be chosen in such a way that the bias and variance ofβ α,L arises from a specific type of structure, potentially decreasing bias without increasing complexity of the model. As a common example, L = D 2 introduces smooth bias structured by low-frequency trigonometric functions.

Bias and variance and the choice of penalty operator
Begin by observing that the penalized estimateβ α,L in (3) is a linear transformation of any solution to the normal equations. Indeed, define X # ≡ X # α,L = (X ′ X +αL ′ L) −1 X ′ and note that ifβ denotes any solution to X ′ Xβ = X ′ y, thenβ α,L = X # Xβ +X # ǫ. The resolution operator X # X reflects the extent to which the estimate in (7) is linearly transformed relative to an exact solution. In particular, Hence bias can be controlled by the choice of L, with an estimate being unbiased whenever Lβ = 0. There is a tradeoff, of course, and equation (10) below quantifies the effect on the variance as determined by W ø (i.e., {w k } n k=n−d+1 ) if Null(L) is chosen to be too large.
More generally, the decompositions in (4) lead to an expression for the resolution matrix as (4), the first p − n rows of W −1 are not used. For notational convenience, definẽ W := W ′−1 (note,W plays a role analogous to V ≡ V ′−1 in the SVD). As before, letW n denote the p × n matrix consisting of the last n columns ofW and note that in equation (4) and so the bias ofβ α,L can be expressed as wherew k is the kth column ofW . In particular, the bias vector (9) is contained in span{w k : µ k = 0}, whereas the estimateβ α,L is in span{w k : σ k = 0}. In the special case that X ′ X is invertible, then [28]). A counterpart is an expression for the variance in terms of the GSVD. Let Σ denote the covariance for ǫ. Then var(β α,L ) = var( An interesting perspective of the bias-variance tradeoff is provided by the relationship between the GS-values in (5) and their role in equations (9) and (10). Moreover, these lead to an explicit expression for the mean squared error (MSE) of a PEER estimate. Since E(β α,L ) = X # Xβ, The GS-vectors {w k } are not necessarily orthogonal, although they are X ′ X-orthogonal; see (6). Consequently, a bound, rather than equality, for the MSE in terms of the GS values/vectors is the best one can do in general: (11)).
As a final remark, recall that one perspective on ridge estimation defines fictitious data from an orthogonal "experiment," represented by an L, and expresses I as I = L ′ L [29]. Regardless of orthogonality this applies to any penalized estimate and L may similarly be viewed as augmenting the data, influencing the estimation process through its eigenstructure; the response, y, is set to zero for these supplementary "data". In this view, equation (3) can be written as Zβ = y where Z = X √ αL and y = y 0 . This formulation proves useful in section 5.3 when assuring that the estimation process is stable with respect to perturbations in X and the choice of L.

Structured penalties
A structured penalty refers to a second term in (1) that involves an operator chosen to encourage certain functional properties in the estimate. A prototypical example is a derivative operator which imposes smoothness via its eigenstructure. Here we describe several examples of structured penalties, including two that were motivated heuristically and implemented without regard to the spectral properties that define their performance. Sections 3.2 and 5.3 provide a complete formulation of their properties as revealed by the GSVD.

The penalty of C. Goutis
The concept of using a penalty operator whose eigenstructure is targeted toward specific properties in the predictors appears implicitly in the work of C. Goutis [14]. This method aimed to account for the "functional nature of the predictors" without oversmoothing and, in essence, considered the inverse of a smoothing penalty. Specifically, if ∆ denotes a discretized second-derivative operator (with some specified boundary conditions), the minimization in (1) was replaced by min Here, the term X∆ ′ ∆β can be viewed as the product of X∆ ′ (derivatives of the predictor curves) and ∆β (derivative of β). Defining γ := ∆ ′ ∆β and seeking a penalized estimate of γ leads toγ In [14], the properties ofγ were conjectured to result from the eigenproperties of (∆ ′ ∆) † . This was explored by ignoring X and plotting some eigenvectors of (∆ ′ ∆) † . The properties of this method become transparent, however, when formulated in terms of the GSVD. That is, let L := ((∆ ′ ∆) † ) 1/2 and note the functions that defineγ are influenced most by the highly oscillatory eigenvectors of L which correspond to its smallest eigenvalues; see equations (5) and (7). This approach was applied in [14] only for prediction and has drawbacks in producing an interpretable estimate, especially for non-smooth predictor curves. The general insight is valid, however, and modifications of this penalty can be used to produce more stable results. The operator (∆ ′ ∆) † essentially reverses the frequency properties of the eigenvectors of ∆ and is an extreme alternative to this smoothing penalty. An eigenanalysis of the pair (X, L), however, suggests penalties that may be more suited to the problem. This is illustrated in Section 7.

Targeted penalties
Given some knowledge about the relevant structure, a penalty can be defined in terms of a subspace containing this structure. For example, suppose β ∈ Q := span{q j } d j=1 in L 2 (I). Set Q = d j=1 q j ⊗ q j and consider the orthogonal projection P Q = QQ † . (Here, q ⊗ q denotes the rank-one operator f → f, q q, for f ∈ L 2 (I).) For L = I − P Q , then β ∈ Null(L) andβ α,L is unbiased. The problem may still be underdetermined so, more pragmatically, define a decomposition-based penalty for some a, b ≥ 0. Heuristically, when a > b > 0 the effect is to move the estimate towards Q by preferentially penalizing components orthogonal to Q; i.e., assign a prior favoring structure contained in the subspace Q. To implement the tradeoff between the two subspaces, we view a and b as inversely related, ab = const. The analytical properties of estimates that arise from this are developed in the next section and illustrated numerically in Section 7. For example, bias is substantially reduced when β ⊂ Q, and equation (19) quantifies the tradeoff with respect to variance when the prior Q is chosen poorly. More generally, one may penalize each subspace differently by defining L = α 1 (I − P Q )L 1 (I − P Q ) + α 2 P Q L 2 P Q , for some operators L 1 and L 2 . This idea could be carried further: for any orthogonal decomposition of L 2 (I) by subspaces Q 1 , . . . , Q J , let P j be the projection onto Q j . Then the multi-space penalty L = J j=1 α j P j leads to the estimatẽ This concept was applied in the context of image recovery (where X represents a linear distortion model for a degraded image y) by Belge et al. [1].
The examples here illustrate ways in which assumptions about the structure of a coefficient function can be incorporated directly into the estimation process. In general, any estimation of β imposes assumptions about its structure (either implicitly or explicitly) and section 3.2 shows that the bias-variance tradeoff involves a choice on the type of bias (spatial structure) as well as the extent of bias (regularization parameter(s)).

Some analytical properties
Any direct comparison between estimates using different penalty operators is confounded by the fact there is no simple connection between the generalized singular values/vectors and the ordinary singular values/vectors. Therefore, we first consider the case of targeted or projection-based penalties (14). Within this class, we introduce a parameterized family of estimates that are comprised of ordinary singular values/vectors. Since the ridge and PCR estimates are contained in (or a limit of) this family, a comparison with some targeted PEER estimates is possible. For more general penalized estimates, properties of perturbations provide some less precise relationships; see proposition 5.6.

Transformation to standard form
We have reason to consider decomposition-based penalties (14) in which L is invertible. In this case, an expression for the estimate does not involve the second term in (8), and decomposing the first term into two parts will be useful. For this, we find it convenient to use the standard-form transformation due to Elden [11] in which the penalty L is absorbed into X. This transformation also provides a computational alternative to the GSVD which, for projection-based penalties in particular, can be less computationally expensive; see, e.g., [25]. By this transformation of X, a general PEER estimate (L = I) can be expressed via a ridge-regression process.
Define the X-weighted generalized inverse of L and the corresponding transformed X as: see [11,19]. In terms of the GSVD components (4), the transformed X is = U ΓV ′ . In particular, the diagonal elements of Γ = SM † are the ordinary singular values of , but in reversed order. Moreover, a PEER estimate can be obtained from a ridge-like penalization process with respect to . That is, for¬ Note that the transformed estimate as given in terms of the GSVD factors is: The regularization parameter, previously denoted by α, can be absorbed into the values a and b, so we will denote this PEER estimateβ α,L simply asβ a,b .
Remark 5.1. When a = b = √ α, this is simply a ridge estimate:β a,b =β α,I . Therefore, the best performance among this family of estimates is as least as good as the performance of ridge, regardless of the choice of Q.

SVD targeted penalties
Consider the special case in which Q is the span of the d largest right singular vectors of an n × p matrix X of rank n. Let X = U 0 D V ′ be an ordinary singular value decomposition where D is a diagonal matrix of singular values. For consistency with the GSVD notation, these will be ordered as 0 ≤ σ 1 ≤ · · · ≤ σ n . As before, the first p − n columns of V are not used. Rather than introduce extra notation, we write X = U DV ′ , letting V denote the n×p whose columns correspond to the singular vectors in D. So now, the last d columns of V correspond to the d largest singular values of X (i.e., Q = V ø ).
We are interested interested in the penalty L = a(I − P Q ) + bP Q , where d = dim(Null(I − P Q )). Similar to before, set o = n − d and define o × o and d × d submatrices, D o and D ø , of D as Here, This decomposition implies that the ridge estimate in (15) is of the following form: setting G = DΛ −1 , denoting its diagonal entries by {γ k }, and defining F = diag{γ 2 k /(γ 2 k + 1)} gives¬ = V F G † U ′ y. Now, By the decomposition (16), This shows that the estimate decomposes as follows.
Theorem 5.2. Let Q be the span of the largest d right singular vectors of X. Set L = a(I − P Q ) + bP Q . Then, in terms of the notation above, the estimateβ a,b decomposes as where the left and right terms are independent of b and a, respectively.
Similar arguments can be used to decompose an estimate for arbitrary Q: In this case, however, all terms are dependent on both a and b. Indeed, using notation as in (8) However, L −1 V = W M † , and the non-orthogonal terms provided by W do not decompose the estimate into terms from the orthogonal sum Q ⊕ Q ⊥ .
The following corollary, along with Remark 5.1, records the manner in which (17) is a family of penalized estimates, parameterized by a, b > 0 and d ∈ {1, ..., n}, that extends some standard estimates. 2. when a > 0 and b = 0,β a,0 is given by (8), which is a sum of PCR and ridge estimates on Q and Q ⊥ , respectively; 3. for each d, the PCR estimateβ d PCR is the limit ofβ a,0 as a → ∞.
In item 2, this estimate is similar to PCR except that a ridge penalty is placed on the leastdominant singular vectors. Under the assumptions here, w k ≡ v k are the ordinary singular vectors of X and the ordinary singular values appear as γ k = σ k /µ k , for µ k > 0. In the second term of (8), the singular vectors are in the null space of L (since b = 0), and so µ k = 0 and σ k = 1, for k = n − d + 1, ..., p. Regarding item 3, although a PCR estimate is not obtained from equation (3) for any L, it is a limit of such estimates.
Other decompositions may be obtained simply by using a permutation, such as Q = ΠV , for some n × n permutation matrix Π. Stein's estimate,β α,S , also fits into this framework as follows.
When X ′ X is nonsingular, thenβ α,S = (X ′ X +αX ′ X) −1 X ′ y (see, e.g., the class 'STEIN' in [10]), and X ′ X = V D ′ DV ′ . Hence this estimate arises from the penalty L S = DV ′ . This is a re-weighted version of L = a(I −P Q ) where d = n, Q = V and the parameter a is replaced by the matrix D. The result is a constant filter factor F = diag{1/(1+ α)}. Using d < n and Q = V d is a natural extension of this idea. More generally, Q may be enriched with functions that span a wider range of structure potentially relevant to the estimate. This concept is illustrated in Section 7.3 where instead of V d , we use a d-dimensional set of experimentally-derived "template" spectra supplemented with their derivatives to define Q.
As an aside, we note that in a different approach to regularization one can define a general family of estimates arising from the SVD by way ofβ h, h 2 )}, and ϕ : Ê + → Ê is an arbitrary continuous function [33]. A ridge estimate is obtained for ϕ(t) = 1/(1+t), and PCR obtained for ϕ(t) = 1/t if t > 1, ϕ(t) = 0 if t ≤ 1 (an L 2 -limit of continuous functions). This is similar to item 3 in Corollary 5.3, but the family of estimatesβ h,ϕ is formulated in terms of functional filter factors rather than explicit penalty operators. Related to this is the fact that the optimal (with respect to MSE) estimate using SVD filter factors is, in the case C = σ ǫ I, expressed }; see the "ideal filter" of [34]. In fact, it's easy to check that this optimal estimate can be obtained asβ OH =β α,L for some L = I. Sincẽ β OH involves knowledge of β, it is not directly obtainable but it points to the optimality of a PEER estimate.

The MSE of some penalized estimates.
Theorem 5.2 is used here to show thatβ a,b can have smaller MSE than the ridge or PCR estimates for a wide range of values of a and/or b. The MSE is potentially decreased further when L is defined by a more general Q. In that case, a general statement is difficult to formulate but Proposition 5.6 confirms that any improvement in MSE is robust to perturbations in L (e.g., general Q) and errors in x.
An immediate consequence of Theorem 5.2 is that the mean squared error for an estimate in this family (17) decomposes into easily-identifiable terms for the bias and variance: The influence of b = 0 on the estimate is now clear: when the numerical rank of X is small relative to d, the σ k 's in the third term decrease and the contribution to the variance from this term increases-the estimate fails for the same reason that ordinary least-squares fails. Any nonzero b stabilizes the estimate in the same way that a nonzero α stabilizes a standard ridge estimate; the decomposition (17) merely re-focuses the penalty. This is illustrated in Section 7 ( Table 1) and in the Appendix (Table 4). Although there are three parameters to consider, the MSE ofβ a,b is relatively insensitive to b > 0 for sufficiently large d. This could be optimized (similar to efforts to optimize the number of principal components) but here we assume some knowledge regarding Q, hence d. Relationships between ridge, PCR and PEER estimates in this family {β a,b } a,b>0 can be quantified more specifically as follows.
Proposition 5.4. Suppose β ∈ Q and fix α > 0. Then for any a > √ α, the ridge estimate satisfies Proof. This follows from the fact that V o ′ β = 0 and so the second term in (19) is zero. Therefore, the contribution to the MSE by the first term is decreased whenever a > √ α.
If β is exactly a sum of the d dominant right singular vectors, A PCR estimate using d terms may perform well, but is not optimal: Proposition 5.5. If β ∈ Q, a sufficient condition for the PCR estimate to satisfy Note that the left side of (20) increases without bound as σ k → 0. Since ||V ø ′ β|| 2 = n k=n−d+1 (v ′ k β) 2 , and since the premise of PCR is that v ′ k β decreases with decreasing σ k , this sufficient condition is entirely plausible.
Proof. The MSE ofβ d PCR consists of the second and third terms of (19): In particular, a sufficient condition for this to exceed MSE(β ∞,b ) is for the variance term to exceed the third and fourth terms of (19):

Its easy to check that this is satisfied when (20) holds.
A comment by Bingham and Larntz [3] on the intensive simulation study of ridge regression in [10] notes that "it is not at all clear that ridge methods offer a clear-cut improvement over [ordinary] least squares except for particular orientations of β relative to the eigenvectors of X ′ X." Equation (19) repeats this observation relating these two classical methods as well as the minor extensions contained in (17). If, on the other hand, "the orientation of β relative to the [v k 's]" is not favorable, i.e., if β is nowhere near the range of V , then a PEER estimate as in (18) is more desirable than (17) (assuming sufficient prior knowledge).
In summary, the family of estimates {β a,b } a,b>0 in (17) represents a hybrid of ridge and PCR estimation. This family-based on the ordinary singular vectors of X-is introduced here to provide a framework within which these two familiar estimates can be compared to (slightly) more general PEER estimates. Direct analytical comparison between general PEER estimates is more difficult since there's no simple relationship between the generalized singular vectors for two different L (including L = I versus L = I). However, it is important that the estimation process be stable with respect to changes in L and/or X. I.e., in going from an estimate in (17) to one in (18), the performance of the estimate should be predictably altered. Given an estimate in Proposition 5.4, if Q is modified and/or X is observed with error, the MSE of the corresponding estimate,β E α,L , should be controlled: for sufficiently small perturbation E, the corresponding estimate MSE(β E α,L ) should be close to MSE(β α,I ). This "stability" is true in general. To see this, recall Z = X ′ √ αL ′ ′ (of rank p) and y = y ′ 0 ′ . Then another way to represent the estimate (3) isβ α,L = Z † y. Let ′ ′ for some n × p and m × p matrices E 1 and E 2 . Set Z E = Z + E and denote the perturbed estimate byβ E α,L = Z † E y. By continuity of the generalized inverse (e.g., [4], Section 1.4), lim ||E||→0 Z † E = Z † if and only if lim ||E||→0 rank(Z E ) = rank(Z). Therefore, provided the rank of Z is not changed by E, and hence MSE(β E α,L ) → MSE(β α,L ) as ||E|| → 0. A more specific bound on the difference of estimates can be obtained under the condition ||Z † ||||E|| < 1 which implies that ||Z † E || < ||Z † || 1−||Z † ||||E|| . This can be used to obtain the following bound.

Tuning parameter selection
Despite our focus on the GSVD, the computation of a PEER estimate in (1) does not, of course, require that this decomposition be computed. Rather, the role of the GSVD has been to provide analytical insight into the role a penalty operator plays in the estimation process. For computation, on the other hand, we have chosen to use a method in which the tuning parameter, α, is estimated as part of the coefficient-function estimation process. Because the choice of tuning parameter is so important, many selection criteria have been proposed, including generalized cross-validation (GCV) [9], AIC and its finite sample corrections [55]. As an alternative to GCV and AIC, a recently-proven equivalence between the penalized least squares estimation and a linear mixed model (LMM) representation [6] can be used. In particular, the best linear unbiased predictor (BLUP) of the response y is composed of the best linear unbiased estimator of the fixed effects and BLUP of the random effects for the given values of the random component variances (see [47] and [6]). Within the LMM framework, restricted maximum likelihood (REML) can be used to estimate the variance components and thus the choice of the tuning parameter, α, which is equal to the ratio of the error variance and the random effects variance [42]. REML-based estimation of the tuning parameter has been shown to perform at least as well as the other criteria and, under certain conditions, it seen to be less variable than GCV-based estimation [41]. In our case, the penalized least-squares criterion (1) is equivalent tõ where β = [β ′ unp β ′ pen ] ′ , the X unp corresponds to the unpenalized part of the design matrix, and X pen to the penalized part.
For simplicity of presentation, we describe the transformation with an invertible L. However, a generalized inverse can be used in case L is not of full rank; see equation (15). Also, to facilitate a straightforward use of existing linear mixed model routines in widely available software packages (e.g., R [37] or SAS software [43]), we transform the coefficient vector β using the inverse of the matrix L. Let X ⋆ = XL and β ⋆ = L −1 β. Then equation (21) can be modified as follows This REML-based estimation of tuning parameters is used in the application of Section 7.3. For estimation of the parameters a, b and α involved in the decomposition-based penalty of equation (17), we view a and b as weights in a tradeoff between the subspaces and assume ab = const. For implementation, we fix one, estimate the other using a grid search, and use REML to estimate α.

Numerical examples
To illustrate algebraic properties given in Section 5, we consider PEER estimation alongside some familiar methods in several numerical examples. Section 7.1 elaborates on the simple example in Section 2.1. These mass spectrometry-like predictors are mathematically synthesized in a manner similar to the study of Reiss and Ogden [40] (see also a numerical study in [48]). Here, β is also synthesized to represent a spectrum, or specific set of bumps. In contrast, Section 7.3 presents a real application to Raman spectroscopy data in which a set of spectra {x i } and nanoparticle concentrations {y i } are obtained from sets of laboratory mixtures. This laboratory-based application is preceded in section 7.2 by a simulation that uses these same Raman spectra. In both Raman examples, targeted penalties (14) are defined using discretized functions q j chosen to span specific subspaces, Q = span{q j } d j=1 . As before, let Q = col[q 1 , ..., q d ] and P Q = QQ † . Each section displays the results from several methods, including derivative-based penalties. Implementing these requires a choice of discretization scheme and boundary conditions which define the operator. We use D 2 where D = [d i,j ] is a square matrix with entries d i,i = 1, d i,i+1 = −1 and d i,j = 0 otherwise. In addition to some standard estimates, sections 7.3 and 7.2 also consider FPCR R , a functional PCR estimate described in [40]. This approach extends the penalized Bspline estimates of [8] and assumes β = Bη where B is an p × K matrix whose columns consist of K B-spline functions and η is a vector of B-spline coefficients. The estimation process takes place in the coefficient space using the penalty L = D 2 applied to η. The FPCR R estimate further assumes β = BV d η (V d as defined in section 2).
Estimation error is defined as mean squared error (MSE) ||β −β α,L || 2 , and the prediction error defined similarly as i |y i −ỹ i | 2 , whereỹ i = x i ,β . Each simulation incorporates response random errors, ǫ i ∼ N (0, σ 2 ǫ ), added to the ith true response, y true i = x i , β . Letting S 2 Y denote the sample variance in the set {y true i } n i=1 , the response random errors, ǫ i , are chosen such that (the squared multiple correlation coefficient of the true model) takes values 0.6 and 0.8. In sections 7.1 and 7.2, tuning parameters are chosen by a grid search. In section 7.3, tuning parameters are chosen using REML, as described in section 6.

Bumps Simulation
Here we elaborate on the simple example of section 2.1. This simulation involves bumpy predictor curves x i (t) with a response y i that depends on the amplitudes x i (t) at some of the bump locations, t = c k , via the regression function β. In particular, where J X = {2, 6, 10, 14, 20, 26, 30}, J β = {6, 14, 26}, a ⋆ are the magnitudes, b ⋆ are the spreads, and c ⋆ are the locations of the bumps. In the first simulation, we set b j = 10000 and c j = 0.004(8j − 1), the same for each curve x i . This mimics, for instance, curves seen in mass spectrometry data.
The assumption J β ⊂ J X simulates a setting in which the response is associated with a subset of metabolite or protein features in a collection of spectra. The a ij 's are from a uniform distribution, and a j = 3, 5, 2 for j = 6, 14, 24, respectively. We consider discretized curves, x i (t), evaluated at p = 250 points, t j , j = 1, . . . , p. The sample size is fixed at n = 50 in each case.
Penalties. We consider a variety of estimation procedures: ridge (L = I), second-derivative (D 2 ), a more general derivative operator (D 2 + a I) and PCR. We also define two decompositionbased penalties (14) formed by specific subspaces Q = span{q j } j∈J for q j of the form q j (t) = a j exp[b j (t − c j )], with c j at all locations seen in the predictors, J V = {2, 6, 10, 14, 20, 26, 30}, or at uniformly-spaced locations, J U = {2, 4, . . . , 30}; denote these penalties by L V and L U , respectively.
Simulation results. The simulation incorporates two sources of noise: (i) response random errors, ǫ i ∼ N (0, σ 2 ǫ ), added to the ith true response so that R 2 = 0.6, 0.8; (ii) measurement error, e i ∼ N p (0, σ 2 e I), added to the ith predictor, x i . To define a signal-to-noise ratio, S/N , set S 2 i := ||x i − µ i || 2 /(p − 1), where µ i is the mean value of x i , and set S 2 X := 1/n i S 2 i . The e i are chosen so that S/N := S X /σ e = 2, 5, 10. Figure 1 shows a few partial sums of (7) for estimates arising from three penalties: D 2 , L = I and L V , when R 2 = 0.8 and S/N = 2. Table 1 gives a summary of estimation errors. The penalty L V , exploiting known structure, performs well in terms of estimation error. Not surprisingly, a penalty that encourages low-frequency singular vectors, D 2 , is a poor choice although D 2 + a I easily improves on D 2 since the GSVs are more compatible with the relevant structure. PCR performs well with estimation errors that can be several times smaller than those of ridge. The number of terms used in PCR ranges here from 8 (S/N = 10) to 25 (S/N = 2). Predictably, PCR performance degrades with decreasing S/N , a property that is less pronounced, or not shared, by other estimates. Performances of L V and L U illustrate properties described in Section 5.3. As S/N → 0, the ordinary singular vectors of X (on which ridge and PCR rely) decreasingly represent the structure in β. The GS vectors of (X, L V ) and (X, L U ), however, retain structure relevant for representing β.   Table 2 summarizes prediction errors. When S/N is large, performance of PCR is comparable with L V and L U , but degrades for low S/N . Here, even D 2 + a I provides smaller prediction errors, in most cases, than ridge, D 2 or PCR. This illustrates the GS vectors role in (13) and reiterates observations in [14].

Raman simulation
We consider Raman spectroscopy curves which represent a vibrational response of laser-excited coorganic/inorganic nanoparticles (COINs). Each COIN has a unique signature spectrum and serves as a sensitive nanotag for immunoassays; see [27,44]. Each spectrum consists of absorbance values measured at p = 600 wavenumbers. By the Beer-Lambert law, light absorbance increases linearly with a COIN's concentration and so a spectrum from a mixture of COINs is reasonably modeled by a linear combination of pure COIN spectra. The data here come from experiments that were  Figure 2: Nine pure COIN spectra, P 1 , ..., P 9 , and a coefficient function, β (each shifted for display). β arises as a solution to the fLM in which y denotes concentrations of P 9 in an in silico mixture of 50 COIN spectra, x i (light gray). This β is used in the simulation study of Section 7.2.
designed to establish the ability of these COINs to measure the existence and abundance of antigens in single-cell assays. Let P 1 , ..., P 10 denote spectra from nine pure COINs and one "blank" (no biochemical material), each normalized to norm one. We form in-silico mixtures as follows: x i = 10 k=1 c i,k P k , i = 1, ..., n, with coefficients {c i,k } generated from a uniform distribution. Figure 2 shows representative spectra from all nine COINs superimposed on a collection of mixture spectra, {x i } 50 i=1 . Included in Figure 2 is the β (dashed curve) used to defined the simulation: y i = x i , β + ǫ, ǫ ∼ N (0, σ 2 ).
In this simulation, we have created a coefficient function which, instead of being modeled mathematically, is a curve that exhibits structure of the type found in Raman spectra. Details on the construction of this β are in Appendix 9.1 so here we simply note that it arises as a ridge estimate from a set of in-silico mixtures of Raman spectra in which one COIN, P 9 , is varied prominently relative to the others. See Figure 2. Motivation for defining β in this way is based on a view that it seems implausible for us to predict the structure of realistic signal in these data and recreate it using polynomials, Gaussians or other analytic functions.
Regardless of its construction, β defines signal that allows us to compute estimation and prediction error. The performances of five methods are summarized in Table 3. Note that although β was constructed as a ridge estimate (using a different set of in-silico mixtures; see Appendix 9.1), the ridge penalty is not necessarily optimal for recovering β. This is because the strictly empirical eigenvectors associated with the new spectra may contain structure not informative regarding y. Also, in these data, the performance of FPCR R is adversely affected by a tendency for the estimate to be smooth; cf., Figure 3. The PEER penalty used here is defined by a decomposition-based operator (17) in which Q is spanned by a 10-dimensional set of pure-COIN spectra (including a blank). The success of such an estimate obviously depends on an informed formation of Q, but as long as the parameter-selection procedure allows for a = b, then the set of possible estimates includes ridge as well as estimates with potentially lower MSE than ridge; see Proposition 5.4.
We note that this simulation may be viewed as inherently unfair since the PEER estimate uses

Raman Application
We now consider spectra representing true antibody-conjugated COINs from nine laboratory mixtures. These mixtures contain various concentrations of eight COINs (of the nine shown in Figure 2). Spectra from four technical replicates in each mixture are included to create a set of n = 36 spectra {x i } n i=1 . We designate P 1 as the COIN whose concentration within each mixture defines y. Assuming a linear relationship between the spectra, {x i }, and the P 1 -concentrations, {y i }, we estimate P 1 . More precisely, we estimate the structure in P 1 that correlates most with its concentrations, as manifest in this set of mixtures. The fLM is a simplistic model of this relationship between the concentration of P 1 and its functional structure, but the physics of this technology imply it is a reasonable starting point.
We present the results of three estimation methods: ridge, FPCR R and PEER. In constructing a PEER penalty, we note that the informative structure in Raman spectra is not that of low-frequency or other easily modeled features, but it may be obtainable experimentally. Therefore, we define L as in (14) in which Q contains the span of COIN template spectra: Q 1 = span{P k } 8 k=1 . However, since a single set of templates may not faithfully represent signal in subsequent experiments (with new measurement errors, background and baseline noise etc), we enlarge Q by adding additional structure related to these templates. For this, set Q 2 = span{P ′ k , P ′′ k } 8 k=1 , where P ′ k denotes the derivative of spectrum P k . (Note, to form Q 2 , scale-based approximations to these derivatives are used since raw differencing of non-smooth spectra introduces noise.) Then set Q = span{Q 1 ∪ Q 2 } and define L = a(I − P Q ) + bP Q .
The regularization parameters in the PEER and ridge estimation processes were chosen using REML, as described in Section 6. For the FPCR R estimate, we used the R-package refund [39] as implemented in [40].
Since β is not known (the model y = Xβ + ǫ is only approximate), we cannot report MSEs for these three methods. However, the structure of P 1 is qualitatively known and by experimental design, y is directly associated with P 1 . The goal here is that of extracting structure of the constituent spectral components as manifest in a linear model. This application is similar to the classic problem of multivariate calibration [5,31] which essentially leads to a regression model using an experimentally-designed set of spectra from laboratory mixes.
The structure in the estimate here is expected to reflect the structure in P 1 that is correlated with P 1 's concentrations, y. The estimate is not, however, expected to precisely reconstruct P 1 since P 1 shares structure with the other COIN spectra not associated with y. See Figure 2 where P 1 is plotted alongside the other COIN spectra. Now, Figure 3 shows plots of the PEER, FPCR R and ridge estimates of the fLM coefficient function. The PEER estimate,β Q , provides an interpretable compromise between ridge, which involves no smoothing, and FPCR R , which appears to oversmooth. For reference, the P 1 spectrum is also plotted along with a mean-adjusted version ofβ Q ,β Q + µ  Finally, we consider prediction for these methods by forming a new set of spectra from different mixture compositions (different concentrations of each COIN) and, additionally, taken from different batches. This "test" set consists of spectra from four technical replicates in each of 15 mixtures forming a set of n = 60 spectra, {x test i } n i=1 . As before, P 1 is the COIN whose concentration within each mixture defines the values {y test i } n i=1 . For the estimates from each of the three methods (shown in Figure 3) we compute the prediction error: The errors for PEER, ridge, and FPCR R estimates are 0.770, 0.752, 2.139, respectively. The ridge estimate here illustrates how low prediction error is not necessarily accompanied by interpretable structure in the estimate (or low MSE) [7].

Discussion
As high-dimensional regression problems become more common, methods that exploit a priori information are increasingly popular. In this regard, many approaches to penalized regression are now founded on the idea of "structured" penalties which impose constraints based on prior knowledge about the problem's scientific setting. There are many ways in which such constraints may be imposed, and we have focused on the algebraic aspects of a penalization process that imposes spatial structure directly into a regularized estimation. This approach fits into the classic framework of L 2 -penalized regression but with an emphasis on the algebraic role that a penalty operator plays to impart structure on the estimate.
The interplay between a structured regularization term and the coefficient-function estimate may not be well understood in part because it is not typically viewed in terms of the generalized singular vectors/values, which is fundamental to this investigation. In particular, any penalized estimate of the form (1) with L = I is intrinsically based on GSVD factors in the same way that many common regression methods (such as PCR, ridge, James-Stein, or partial least squares) are intrinsically based on SVD factors. Just as the basics of the ubiquitous SVD are important to understanding these methods, we have aspired to established the basics of the GSVD as it applies to a this general penalized regression setting and to illustrate how the theory underlying this approach can be used inform the choice of penalty operator.
Toward this goal the presentation emphasizes the transparency provided by the partially-empirical eigenvector expansion (7). Properties of the estimate's variance and bias are determined explicitly by the generalized singular vectors whose structure is determined by the penalty operator. We have restricted attention to additive constraints defined by penalty operators on L 2 in order to retain the direct algebraic connection between the eigenproperties of the operator pair (X, L) and the spatial structure ofβ α,L . Intuitively, the structure of the penalty's least-dominant singular vectors should be commensurate with the informative structure of β. The actual effect a penalty has on the properties of the estimate can be quantified in terms of the GSVD vectors/values. This perspective differs from popular two-stage signal regression methods in which estimation is either preceded by fitting the predictors to a set of (external) basis functions or is followed by a step that smooths the estimate [8,21,30,38,40]. Instead, structure (smoothness or otherwise) is imposed directly into the estimation process. The implementation of a penalty that incorporates structure less generic than smoothness (or sparseness) requires some qualitative knowledge about spatial structure that is informative. Clearly this is not possible in all situations, but our presentation has focused on how a functional linear model may provide a rigorous and analytically tractable way to take advantage of such knowledge when it exists. 9 Appendix 9.1 Defining β for the simulation in section 7.2 This simulation is motivated by an interest in constructing a plausibly realistic β whose structure is naturally derived by the scientific setting involving Raman signatures of nanoparticles. Although one could model a β mathematically using, say, polynomials or Gaussian bumps (cf., Appendix A.2), such a simulation would be detached from the physical nature of this problem. Instead, we construct a coefficient function that genuinely comes from a functional linear model with Raman spectra as predictors.
We first generate in-silico mixtures of COIN spectra as x o i = 9 k=1 c i,k P k , i = 1, ..., 50, where c i,k ∼ unif[0, 1]. Designating P 9 as the COIN of interest, we define response values that correspond to the "concentration" of P 9 by setting y o i := 3 c i,9 , i = 1, ..., n. The factor of 3 imposes a strong association between P 9 and the response. Now, the example in section 7.2 aims to estimate a coefficient function that truly comes from a solution to a linear model. However, the equation y o = X o β has infinitely many solutions (where X o is the matrix whose ith row is x o i ), so we must we must regularize the problem to obtain a specific β. For this, we simply use a ridge penalty and designate the resulting solution to be β. This is shown by the dotted curve in Figure 2 and is qualitatively similar to P 9 .
We note that the simulation in section 7.2 uses the same set of COINs, but a new set of in-silico mixture spectra (i.e., a new set of {c i,k } ∼ Unif[0, 1]). In addition, a small amount of measurement error was added, as in section 7.1, to each spectrum during the simulation.

Frequency domain simulation
We display results from a study that mimics the scenario of simulations studied by Hall and Horowitz [16]. We illustrate, in particular, properties of the MSE discussed following equation (19) in section 5.3 relating to b = 0. In fact, we consider the more general scenario in which Q is not constructed from empirical eigenvectors (as in PCR and ridge), but is defined by a prespecified envelope of frequencies.
Penalties. We consider properties of estimates from a variety of penalties: ridge (L = I), D 2 , D 2 + aI, and PCR 1 . In addition, targeted penalties of the form L = I − P Q , are defined by the specified subspaces Q = span{φ j } j∈J , for φ j defined above. Specifically, we use J = J F = {j = 5, ..., 17} (a tight envelope of frequencies) to define L F , and J = J G = {j = 4, ..., 20} (a less focused span of frequencies) to define L G . The operator D 2 + aI simply serves to illustrate the role of higher-frequency singular vectors as discussed in Section 4.1. In the simulations, the coefficient a in D 2 + aI was chosen simultaneously with α via a two-dimensional grid search.
Simulation results. Table 4 summarizes estimation results for all six penalties and two sample sizes, n = 50, 200. The prediction results for these estimates are in Table 5. These are reported for S/N = 10, 5 and R 2 = 0.8, 0.6. The number of terms in the PCR estimate was optimized and ranged from 19 to 25 when R 2 = 0.8 and decreased with decreasing R 2 . Analogously, one could optimize over the dimension of Q (to implement a truncated GSVD), but the purpose here is illustrative while in practice a more robust approach would emply a penalty of the form (14). Errors obtained with ridge and P CR are small, as expected, since the structure of β in this example is consistent with the structure represented in the singular vectors, v k . Therefore, even though the relationship between the y i and x i degrades (indeed, even as R 2 → 0), these estimates are comprised of vectors that generally capture structure in β since it is strongly represented by the dominant eigenstructure of X. The second-derivative penalty, D 2 , produces the worst estimate in each of the scenarios due to oversmoothing. Note D 2 + a I improves on D 2 , yet it is still not optimal for the range of frequencies in β.
Regarding L G , the MSE gets worse as S/N increases. Indeed, here Q is fixed and relatively large and since the σ k decay faster when S/N is big, this leads to rank deficiency and large variance; see equation (19) (note, this only applies approximately since Q does not consist of ordinary SVs). In our previous examples, this is stabilized by a b > 0.
The problems of estimation and prediction have different properties [7]; good prediction may be obtained even with a poor estimate, as seen in Table 5. The estimate from L Da is generally poor relative to others (as measured by the L 2 -norm), but its prediction error is comparable with other methods and is best among the non-targeted penalization methods. This is consistent with the outcome described by C. Goutis [14] where (derivatives of) the predictor curves contain sharp features and so standard least-squares regularization (OLS, PCR, ridge, etc.) perform worse than a PEER estimate which imposes a greater emphasis on "regularly oscillatory but not smooth components"; see section 4.1.