A GENERALIZED WEIGHTED OPTIMIZATION METHOD
Year: 2021

Abstract:
The generalization capacity of various machine learning models exhibits different phenomena in the underand over-parameterized regimes. In this paper, we focus on regression models such as feature regression and kernel regression and analyze a generalized weighted least-squares optimization method for computational learning and inversion with noisy data. The highlight of the proposed framework is that we allow weighting in both the parameter space and the data space. The weighting scheme encodes both a priori knowledge on the object to be learned and a strategy to weight the contribution of different data points in the loss function. Here, we characterize the impact of the weighting scheme on the generalization error of the learning method, where we derive explicit generalization errors for the random Fourier feature model in both the underand over-parameterized regimes. For more general feature maps, error bounds are provided based on the singular values of the feature matrix. We demonstrate that appropriate weighting from prior knowledge can improve the generalization capability of the learned model.

heading: 1 INTRODUCTION
text: Given N data pairs {xj , yj}Nj=1, where xj ∈ R, yj ∈ C, j = 1, . . . , N , we are interested in learning a random Fourier feature (RFF) model (Rahimi & Recht, 2008; Liao et al., 2020; Xie et al., 2020)
fθ(x) = P−1∑ k=0 θke ikx, x ∈ [0, 2π], (1)
where P ∈ N is a given positive integer and we used the short-hand notation θ := (θ0, · · · , θP−1)T with the superscript T denoting the transpose operation.
This exact model as well as its generalization to more complicated setups have been extensively studied; see for instance Liao & Couillet (2018); Shahrampour & Kolouri (2019); d’Ascoli et al. (2020); Li et al. (2020); Özcelikkale (2020); Liu et al. (2020; 2021) and references therein. While this model may seem to be overly simplified from a practical perspective for many real-world applications, it serves as a prototype for theoretical understandings of different phenomena in machine learning models (Sriperumbudur & Szabo, 2015; Belkin et al., 2020; Li et al., 2021a).
A common way to computationally solve this learning problem is to reformulate it as an optimization problem where we find θ by minimizing the model and data mismatch for a given dataset. In this paper, we assume that the training data are collected on a uniform grid of x over the domain [0, 2π]. That is, {xj = 2πjN } N−1 j=0 . Let ωN = exp( 2πi N ) where i is the imaginary unit. We introduce Ψ ∈ CN×P to be the feature matrix with elements (Ψ)jk = (ωN )
jk, 0 ≤ j ≤ N − 1, 0 ≤ k ≤ P − 1 . Based on the form of fθ(x) in (1), we can then write the 2-norm based data mismatch into the form∑N−1 j=0 |fθ(xj) − yj |2 = ‖Ψθ − y‖22 where the column data vector y = (y0, · · · , yN−1)T. The learning problem is therefore recast as a least-squares optimization problem of the form
θ̂ = arg min θ
‖Ψθ − y‖22, (2)
assuming that a minimizer does exist, especially when we restrict θ to an appropriate space.
In a general feature regression problem, the Fourier feature {eikx}P−1k=0 is then replaced with a different feature model {ϕk(x)}P−1k=0 , while the least-squares form (2) remains unchanged except that the entries of the matrix Ψ is now Ψjk = ϕk(xj). We emphasize that this type of generalization will be discussed in Section 5. Moreover, we remark that this least-squares optimization formulation is a classical computational inversion tool in solving the general linear inverse problems of the form Ψθ = y; see for instance Engl et al. (1996); Tarantola (2005) and references therein.
Previous work on weighted optimization for feature and kernel learning. Xie et al. (2020) studied the fitting problem for this model under the assumption that the coefficient vector θ is sampled from a distribution with the property that γ is a positive constant,
Eθ[θ] = 0, Eθ[θθ∗] = cγΛ−2γ[P ] , (3)
where the superscript ∗ denotes the Hermitian transpose and the diagonal matrix Λ[P ] has diagonal elements (Λ[P ])kk = tk = 1 + k, k ≥ 0. That is,
Λ[P ] = diag{t0, t1, t2, . . . , tk, . . . , tP−1}, tk := 1 + k . (4)
The subscript [P ] indicates that Λ[P ] is a diagonal submatrix of Λ that contains its element indexed in the set [P ] := {0, 1, · · · , P − 1}. The normalization constant cγ = 1/( ∑P−1 k=0 (1 + k)
−2γ) is only selected so that Eθ[‖θ‖2] = 1. It does not play a significant role in the rest of the paper. The main assumption in (3) says that statistically, the signal to be recovered has algebraically decaying Fourier coefficients. This is simply saying that the target function we are learning is relatively smooth, which is certainly the case for many functions as physical models in practical applications.
It was shown in Xie et al. (2020) that, to learn a model with p ≤ P features, it is advantageous to use the following weighted least-squares formulation
θ̂p = Λ −β [p] ŵ, with ŵ = arg minθ ‖Ψ[N×p]Λ −β [p] w − y‖ 2 2, (5)
when the learning problem is overparameterized, i.e., p > N . Here, Ψ[N×p] ∈ CN×p is the matrix containing the first p columns of Ψ, and β > 0 is some pre-selected exponent that can be different from the γ in (3). To be more precise, we define the the generalization error of the learning problem
Eβ(P, p,N) := Eθ [ ‖fθ(x)− fθ̂p(x)‖ 2 L2([0,2π]) ] = Eθ [ ‖θ̂p − θ‖22 ] , (6)
where the equality comes from the Parseval’s identity, and θ̂p is understood as the vector (θTp , 0, · · · , 0)T so that θ and θ̂p are of the same length P . The subscript θ in Eθ indicates that the expectation is taken with respect to the distribution of the random variable θ. It was shown in Xie et al. (2020) that the lowest generalization error achieved from the weighted least-squares approach (5) in the overparameterized regime (p > N ) is strictly less than the lowest possible generalization error in the underparameterized regime (p ≤ N ). This, together with the analysis and numerical evidence in previous studies such as those in Belkin et al. (2019; 2020), leads to the understanding that smoother approximations (i.e., solutions that are dominated by lower Fourier modes) give better generalization in learning with the RFF model (1).
Main contributions of this work. In this work, we analyze a generalized version of (5) for general feature regression from noisy data. Following the same notations as before, we introduce the following weighted least-squares formulation for feature regression:
θ̂ δ p = Λ −β [p] ŵ, with ŵ = arg min w ‖Λ−α[N ]
( Ψ[N×p]Λ −β [p] w − y δ ) ‖22, (7)
where the superscript δ on y and θ̂p denotes the fact that the training data contain random noise of level δ (which will be specified later). The exponent α is pre-selected and can be different from β. While sharing similar roles with the weight matrix Λ−β[P ] , the weight matrix Λ −α [N ] provides us the additional ability to deal with noise in the training data. Moreover, as we will see later, the weight matrix Λ−α[N ] does not have to be either diagonal or in the same form as the matrix Λ −β [p] ; the current form is to simplify the calculations for the RFF model. It can be chosen based on the a priori information we have on the operator Ψ as well as the noise distribution of the training data.
The highlight and also one of the main contributions of our work is that we introduce a new weight matrix Λ−α[N ] that emphasizes the data mismatch in terms of its various modes, in addition to Λ −β [p] , the weight matrix imposed on the unknown feature coefficient vector θ. This type of generalization has appeared in different forms in many computational approaches for solving inverse and learning problems where the standard 2-norm (or `2 in the infinite-dimensional setting) is replaced with a weighted norm that is either weaker or stronger than the unweighted 2-norm.
In this paper, we characterize the impact of the new weighted optimization framework (7) on the generalization capability of various feature regression and kernel regression models. The new contributions of this work are threefold. First, we discuss in detail the generalized weighted leastsquares framework (7) in Section 2 and summarize the main results for training with noise-free data in Section 3 for the RFF model in both the overparameterized and the underparameterized regimes. This is the setup considered in Xie et al. (2020), but our analysis is based on the proposed weighted model (7) instead of (5) as in their work. Second, we provide the generalization error in both two regimes for the case of training with noisy data; see Section 4. This setup was not considered in Xie et al. (2020), but we demonstrate here that it is a significant advantage of the weighted optimization when data contains noise since the weighting could effectively minimize the influence of the noise and thus improve the stability of feature regression. Third, we extend the same type of results to more general models in feature regression and kernel regression that are beyond the RFF model, given that the operator Ψ satisfies certain properties. In the general setup presented in Section 5, we derive error bounds in the asymptotic limit when P , N , and p all become very large. Our analysis provides some guidelines on selecting weighting schemes through either the parameter domain weighting or the data domain weighting, or both, to emphasize the features of the unknowns to be learned based on a priori knowledge.

heading: 2 GENERALIZED WEIGHTED LEAST-SQUARES FORMULATION
text: There are four essential elements in the least-squares formulation of the learning problem: (i) the parameter to be learned (θ), (ii) the dataset used in the training process (y), (iii) the feature matrix (Ψ), and (iv) the metric chosen to measure the data mismatch between Ψθ and y.
Element (i) of the problem is determined not only by the data but also by a priori information we have. The information encoded in (3) reveals that the size (i.e., the variance) of the Fourier modes in the RFF model decays as fast as (1 + k)−2γ . Therefore, the low-frequency modes in (1) dominate high-frequency modes, which implies that in the learning process, we should search for the solution vectors that have more low-frequency components than the high-frequency components. The motivation behind introducing the weight matrix Λ−β[p] in (5) is exactly to force the optimization algorithm to focus on admissible solutions that are consistent with the a priori knowledge given in (3), which is to seek θ whose components |θk|2 statistically decay like (1 + k)−2β . When the problem is formally determined (i.e., p = N ), the operator Ψ is invertible, and the training data are noise-free, similar to the weight matrix Λ−β[p] , the weight matrix Λ −α [N ] does not change the solution of the learning problem. However, as we will see later, these two weight matrices do impact the solutions in various ways under the practical setups that we are interested in, for instance, when the problem is over-parameterized or when the training data contain random noise.
The weight matrix Λ−α[N ] is introduced to handle elements (ii)-(iv) of the learning problem. First, since Λ−α[N ] is directly applied to the data y
δ , it allows us to suppress (when α > 0) or promote (when α < 0) high-frequency components in the data during the training process. In particular, when transformed back to the physical space, the weight matrix Λ−α[N ] with α > 0 corresponds to a smoothing convolutional operator whose kernel has Fourier coefficients decaying at the rate k−α. This operator suppresses high-frequency information in the data. Second, Λ−α[N ] is also directly applied to Ψθ. This allows us to precondition the learning problem by making Λ−α[N ]Ψ a betterconditioned operator (in an appropriate sense) than Ψ, for some applications where the feature matrix Ψ has certain undesired properties. Finally, since Λ−α[N ] is applied to the residual Ψθ − y, we can regard the new weighted optimization formulation (7) as the generalization of the classic leastsquares formulation with a new loss function (a weighted norm) measuring the data mismatch.
Weighting optimization schemes such as (7) have been studied, implicitly or explicitly, in different settings (Needell et al., 2014; Byrd & Lipton, 2019; Engquist et al., 2020; Li, 2021; Yang et al., 2021). For instance, if we take β = 0, then we have a case where we rescale the classical leastsquares loss function with the weight Λ−α[N ]. If we take α = 1, then this least-squares functional is equivalent to the loss function based on the H−1 norm, instead of the usual L2 norm, of the mismatch between the target function fθ(x) and the learned model fθ̂(x). Based on the asymptotic equivalence between the quadratic Wasserstein metric and the H−1 semi-norm (on an appropriate functional space), this training problem is asymptotically equivalent to the same training problem based on a quadratic Wasserstein loss function; see for instance Engquist et al. (2020) for more detailed illustration on the connection. In the classical statistical inversion setting, Λ2α plays the role of the covariance matrix of the additive Gaussian random noise in the data (Kaipio & Somersalo, 2005). When the noise is sampled from mean-zero Gaussian distribution with covariance matrix Λ2α, a standard maximum likelihood estimator (MLE) is often constructed as the minimizer of (Ψθ − y)∗Λ−2α[N ] (Ψθ − y) = ‖Λ −α [N ](Ψθ − y)‖ 2 2.
The exact solution to (7), with X+ denoting the Moore–Penrose inverse of operator X, is given by
θ̂ δ p = Λ −β [p] ( Λ−α[N ]Ψ[N×p]Λ −β [p] )+ Λ−α[N ]y δ . (8)
In the rest of this paper, we analyze this training result and highlight the impact of the weight matrices Λ−α[N ] and Λ −β [N ] in different regimes of the learning problem. We reproduce the classical bias-variance trade-off analysis in the weighted optimization framework. For that purpose, we utilize the linearity of the problem to decompose θ̂ δ
p as
θ̂ δ p = Λ −β [p] (Λ −α [N ]Ψ[N×p]Λ −β [p] ) +Λ−α[N ]y + Λ −β [p] (Λ −α [N ]Ψ[N×p]Λ −β [p] ) +Λ−α[N ](y δ − y), (9)
where the first part is simply θ̂p, the result of learning with noise-free data, while the second part is the contribution from the additive noise. We define the generalization error in this case as
Eδα,β(P, p,N) = Eθ,δ [ ‖fθ(x)− fθ̂δp(x)‖ 2 L2([0,2π]) ] = Eθ,δ [ ‖θ̂ δ p − θ̂p + θ̂p − θ‖22 ] , (10)
where the expectation is taken over the joint distribution of θ and the random noise δ. By the standard triangle inequality, this generalization error is bounded by sum of the generalization error from training with noise-free data and the error caused by the noise. We will use this simple observation to bound the generalization errors when no exact formulas can be derived. We also look at the variance of the generalization error with respect to the random noise, which is
Varδ(Eθ[‖θ̂ δ − θ‖22]) := Eδ[(Eθ[‖θ̂ δ − θ‖22]− Eθ,δ[‖θ̂ δ − θ‖22])2] . (11)
In the rest of the work, we consider two parameter regimes of learning: (i) In the overparameterized regime, we have the following setup of the parameters:
N < p ≤ P, and, P = µN, p = νN for some µ, ν ∈ N s.t. µ ≥ ν 1 . (12)
(ii) In the underparameterized regime, we have the following scaling relations:
p ≤ N ≤ P, and, P = µN for some µ ∈ N. (13)
The formally-determined case of p = N ≤ P is included in both the overparameterized and the underparameterized regimes. We make the following assumptions throughout the work:
(A-I) The random noise δ in the training data is additive in the sense that yδ = y + δ. (A-II) The random vectors δ and θ are independent. (A-III) The random noise δ ∼ N (0, σI[P ]) for some constant σ > 0.
While assumptions (A-I) and (A-II) are essential, assumption (A-III) is only needed to simplify the calculations. Most of the results we obtain in this paper can be reproduced straightforwardly for the random noise δ with any well-defined covariance matrix.

heading: 3 GENERALIZATION ERROR FOR TRAINING WITH NOISE-FREE DATA
text: We start with the problem of training with noise-free data. In this case, we utilize tools developed in Belkin et al. (2020) and Xie et al. (2020) to derive exact generalization errors. Our main objective is to compare the difference and similarity of the roles of the weight matrices Λ−α[N ] and Λ −β [p] .
We have the following results on the generalization error. The proof is in Appendix A.1 and A.2. Theorem 3.1 (Training with noise-free data). Let δ = 0, and θ be sampled with the properties in (3). Then the generalization error in the overparameterized regime (12) is:
E0α,β(P, p,N) = 1+cγ p−1∑ j=0 t−2γj −2cγ N−1∑ k=0
∑ν−1 η=0 t
−2β−2γ k+Nη∑ν−1
η=0 t −2β k+Nη
+cγ N−1∑ k=0
( ∑ν−1 η=0 t −4β k+Nη)( ∑µ−1 η=ν t −2γ k+Nη)
( ∑ν−1 η=0 t −2β k+Nη) 2 ,
(14) and the generalization error in the underparameterized regime (13) is:
E0α,β(P, p,N) = cγ P−1∑ j=N t−2γj + cγ p−1∑ k=0 µ−1∑ η=1 t−2γk+Nη− cγ N−1∑ k=p µ−1∑ η=1 t−2γk+Nη +N N−p−1∑ i,j=0 ẽ (N) ij ê (N) ji ΣiiΣjj ,
(15)
where {tj}P−1j=0 and cγ are those introduced in (3) and (4), while
ẽ (N) ij = ∑N−1 k=0 t 2α k UkiUkj , ê (N) ij = ∑N−p−1 k′=0 (cγt −2γ p+k′+χp+k′)V ik′Vjk′ , 0 ≤ i, j ≤ N−p−1 ,
with UΣV∗ being the singular value decomposition of Λα[N ]Ψ[N ]\[p] and {χm}N−1m=0 defined as
χm = ∑N−1 k=0 (∑µ−1 η=1 t −2γ k+Nη )( 1 N ∑N−1 j=0 ω (m−k)j N ) , 0 ≤ m ≤ N − 1 .
We want to emphasize that the generalization errors we obtained in Theorem 3.1 are for the weighted optimization formulation (7) where we have an additional weight matrix Λ−α[N ] compared to the formulation in Xie et al. (2020), even though our results look similar to the previous results of Belkin et al. (2020) and Xie et al. (2020). Moreover, we kept the weight matrix Λ−β[p] in the underparameterized regime, which is different from the setup in Xie et al. (2020) where the same weight matrix was removed in this regime. The reason for keeping Λ−β[p] in the underparameterized regime will become more obvious in the case of training with noisy data, as we will see in the next section.
Here are some key observations from the above results, which, we emphasize again, are obtained in the setting where the optimization problems are solved exactly, and the data involved contain no random noise. The conclusion will differ when data contain random noise or when optimization problems cannot be solved exactly.
First, the the weight matrix Λ−β[p] only matters in the overparameterized regime while Λ −α [N ] only matters in the underparameterized regime. In the overparameterized regime, the weight Λ−β[p] forces the inversion procedure to focus on solutions that are biased toward the low-frequency modes. In the underparameterized regime, the matrix Λ−α[N ] re-weights the frequency content of the “residual”
(data mismatch) before it is backprojected into the learned parameter θ̂. Using the scaling P = µN and p = νN in the overparameterized regime, and the definition of cγ , we can verify that when α = β = 0, the generalization error reduces to
E00,0(P, p,N) = 1− N
p +
2N p cγ ∑P−1 j=p t −2γ j = 1 + N p − 2N p cγ ∑p−1 j=0 t −2γ j . (16)
This is given in (Xie et al., 2020, Theorem 1). We mention, however, even in the case of α = 0, when β 6= 0, our result in (14) offered a minor correction to its equivalence in Xie et al. (2020). Second, when the learning problem is formally determined, i.e., when p = N , neither weight matrices play a role when the training data are generated from the true model with no additional random noise and the minimizer can be found exactly. The generalization error simplifies to
E0α,β(P, p,N) = 2cγ ∑P−1 j=p t −2γ j . (17)
This is not surprising as, in this case, Ψ is invertible (because it is a unitary matrix scaled by the constant N ). The true solution to Ψθ = y is simply θ = Ψ−1y. The weight matrices in the optimization problem are invertible and therefore do not change the true solution of the problem. The generalization error, in this case, is therefore only due to the Fourier modes that are not learned from the training data, i.e., modes p to P − 1. While it is obvious from the formulas for the generalization error that the weight matrix Λ−α[N ] indeed plays a role in the underparameterized regime, we show in Figure 1a the numerical calculation of the singular value decomposition of the matrix Λα[N ]Ψ[N ]\[p] for the case of (N, p) = (1024, 512). The impact of α can be seen by comparing the singular values to their correspondence in the α = 0 case (where all the singular values are the same and equal to √ N ). When the system size is large (in this case N = 1024), even a small α (e.g., α = 0.4) can significantly impact the result. In Figure 1b, we plot the theoretical prediction of E0α,β in Theorem 3.1 to demonstrate the doubledescent phenomenon observed in the literature on the RFF model. We emphasize again that in this particular noise-free setup with perfectly solved minimization problem by the pseudoinverse, Λ−α[N ] only plays a role in the underparameterization regime as can be seen from the double-descent curves.
Selecting p to minimize generalization error. It is clear (and also expected) from (16) and (17) that, in the cases of p = N or α = β = γ = 0, the generalization error decreases monotonically with respect to the number of modes learned in the training process. One should then learn as many Fourier coefficients as possible. When p 6= N or α, β 6= 0, E0α,β(P, p,N), for fixed P and N , does not change monotonically with respect to p anymore. In such a situation, we can choose the p values that minimize the generalization error and perform learning with these p values.

heading: 4 ERROR BOUNDS FOR TRAINING WITH NOISY DATA
text: In this section, we study the more realistic setting of training with noisy data. The common practice is that when we solve the minimization problem in such a case, we should avoid overfitting the model to the data by stopping the optimization algorithm early at an appropriate level of values depending on the noise level in the data, but see Bartlett et al. (2020); Li et al. (2021b) for some analysis in the direction of “benign overfitting”. When training with noisy data, the impact of the weight matrices Λ−α[N ] and Λ −β [p] on the generalization error becomes more obvious. In fact, both weight matrices play non-negligible roles in the overparameterized and the underparameterized regimes, respectively.
We start with the circumstance where we still match the data perfectly for each realization of the noise in the data. The result is summarized as follows.
Lemma 4.1 (Training with noisy data: exact error). Under the assumptions (A-I)-(A-III), the generalization error Eδα,β(P, p,N) is given as
Eδα,β(P, p,N) = E0α,β(P, p,N) + Enoise(P, p,N) ,
where E0α,β(P, p,N) is the generalization error from training with noise-free data given in Theorem 3.1 and Enoise(P, p,N) is the error due to noise. The error due to noise and the variance of the generalization error with respect to noise are respectively Enoise(P, p,N) = σ2 N−1∑ k=0 ∑ν−1 η=0 t −4β k+Nη[∑ν−1
η=0 t −2β k+Nη
]2 , Varδ(Eθ[‖θδ−θ‖22]) = 2σ4N2 N−1∑ k=0
[∑ν−1 η=0 t −4β k+Nη ]2 [∑ν−1
η=0 t −2β k+Nη ]4 . (18)
in the overparameterized regime (12), and
Enoise(P, p,N) = σ2 (2p−N
N + N−p−1∑ j=0 ẽ (N) jj Σ2jj ) , p > N/2,
Varδ ( Eθ[‖θδ − θ‖22] ) = 2σ4 (2p−N N2 + N−p−1∑ i,j=0 ẽ (N) ij ẽ (N) ji Σ2iiΣ 2 jj ) ,
(19)
in the underparameterized regime (13).
The proof of this lemma is documented in Appendix A.3. Note that due to the assumption that the noise δ and the coefficient θ are independent, the noise-averaged generalization errors are split exactly into two separate parts: the part due to θ and the part due to δ. The coupling between them disappears. This also leads to the fact that the variance of the generalization error with respect to noise only depends on the noise distribution (instead of both the noise variance and the θ variance).
For the impact of noise on the generalization error, we see again that the impact of Λ−β[P ] is only seen in the overparameterized regime while that of Λ−α[N ] is only seen in the underparameterized regime. This happens because we assume that we can solve the optimization problem exactly for each realization of the noisy data. In the overparameterized regime, when β = 0, we have that Enoise = Nσ2/p as expected. In such a case, the variance reduces to Varδ(Eθ[‖θδ − θ‖22]) = 2Nσ4/p2. In the underparameterized regime, when α = 0, we have that Enoise = pσ2/N . The corresponding variance reduces to Varδ(Eθ[‖θδ − θ‖22]) = 2pσ4/N2. In the limit when p → N , that is, in the formally determined regime, the mean generalization error due to random noise (resp. the variance of the generalization error with respect to noise) converges to the same value Enoise(P, p,N) = σ2 (resp. Varδ(Eθ[‖θδ − θ‖22]) = 2σ4/N ) from both the overparameterized and the underparameterized regimes. This is the classical result in statistical learning theory (Kaipio & Somersalo, 2005).
The explicit error characterization above is based on the assumption that we solve the optimization exactly by minimizing the mismatch to 0 in the learning process. In practical applications, it is often the case that we stop the minimization process when the value of the loss function reaches a level that is comparable to the size of the noise in the data (normalized by the size of the target function, for instance). We now present the generalization error bounds for such a case. Theorem 4.2 (Training with noisy data: error bounds). In the same setup as Lemma 4.1, under the scaling p = N P →∞, we have the following general bound on the generalization error:
Eδα,β(P, p,N) . p−2α−1Eδ[‖δ‖22,Λ−α [N] ] + p−2β(1− p−2γ+1)Eθ[‖θ‖22,Λ−β [p] ].
When α ≥ −1/2, the error decreases monotonically with respect to p. When α < −1/2, the bound is minimized by selecting
p ∼ (Eδ[‖δ‖22,Λ−α [N] ]−1Eθ[‖θ‖22,Λ−β [p] ]) 1 2(α−β) , (20)
in which case we have
Eδα,β(P, p,N) . Eθ[‖θ‖22,Λ−β [p]
] −2α+1 2(α−β)Eδ[‖δ‖22,Λ−α
[N]
] 2β+1 2(α−β) . (21)
Note that the calculations in (20) and (21) are only to the leading order. We neglected the contribution from the term involving γ. Also, the quantities in expectations can be simplified. We avoided doing so to make these quantities easily recognizable as they are the main quantities of interests.
In the underparameterized regime, that is, the regime of reconstruction, classical statistical inversion theory shows that it is statistically beneficial to introduce the weight matrix Λ−α[N ] that is related to the covariance matrix of the noise in the data (Kaipio & Somersalo, 2005) (see also Bal & Ren (2009) for an adaptive way of adjusting the weight matrix for some specific applications). Our result here is consistent with the classical result as we could see that if we take Λ−α[N ] to be the inverse of the covariance matrix for the noise distribution, the size of the noise contribution in the generalization error is minimized.
Theorem 4.2 shows that, in the case of training with noisy data, the α parameter can be tuned to reduce the generalization error of the learning process just as the β parameter in the weight matrix Λ−β[P ] . Moreover, this result can serve as a guidance on the selection of the number of features to be pursued in the training process to minimize the generalization error, depending on the level of random noise in the training data as well as the regime of the problem.

heading: 5 EXTENSION TO GENERAL FEATURE REGRESSION
text: The results in the previous sections, even though are obtained for the specific form of Fourier feature regression, also hold in more general settings. Let fθ be a general random feature model of the form
fθ(x) = ∑P−1 k=0 θkϕk(x), x ∈ X, (22)
constructed from a family of orthonormal features {ϕk}Pk=0 in L2(X). Let Ψ be the feature matrix constructed from the dataset {xj , yj}N−1j=0 :
(Ψ)jk = ϕk(xj), 0 ≤ j, k ≤ P − 1 .
We can then apply the same weighted optimization framework (7) to this general model. As in the case of RFF model, the data in the learning process allows decomposition
yδ = Ψ[N ]×[p]θ[p] + Ψ[N ]×([P ]\[p])θ[P ]\[p] + δ .
The learning process only aims to learn the first pmodes, i.e., θ[p], indicating that the effective noise that we backpropagate into θ[p] in the learning process is Ψ[N ]×([P ]\[p])θ[P ]\[p] + δ. The frequency contents of this effective noise can come from the true noise δ, the part of Fourier modes that we are not learning (i.e., θ[P ]\[p]), or even the feature matrix Ψ[N ]×([P ]\[p]). For the RFF model, the feature matrix Ψ[N ]×[N ] is unitary after being normalized by 1/ √ N . Therefore, all its singular values are homogeneously √ N . For learning with many other features in practical applications, we could have feature matrices Ψ[N ]×[N ] with fast-decaying singular values. This, on one hand, means that Ψ[N ]×([P ]\[p])θ[P ]\[p] will decay even faster than θ[P ]\[p], making its impact on the learning process smaller. On the other hand, Ψ[N ]×[N ] having fast-decaying singular values will make learning θ[p] harder (since it is less stable).
The impact of the weighting scheme on the generalization error as an expectation over θ and δ (defined in Lemma 4.1) is summarized in Theorem 5.1. Its proof is in Appendix A.4. Note that here we do not assume any specific structure on the distributions of θ and δ except there independence. Theorem 5.1. Let Ψ = UΣV∗ be the singular value decomposition of Ψ. Under the assumptions (A-I)-(A-III), assume further that Σkk ∼ t−ζk for some ζ > 0. Then the generalization error of training, using the weighted optimization scheme (7) with Λ−α[N ] replaced with Λ −α [N ]U
∗ and Λ−β[P ] replaced with VΛ−β[P ] , satisfies
Eδα,β(P, p,N) . p2(ζ−α)Eδ[‖δ‖22,Λ−α [N] ] + p−2βEθ[‖θ‖22,Λ−β [p] ] ,
in the asymptotic limit p = N = P →∞. The bound is minimized, as a function of p, when
p ∼ (Eδ[‖δ‖22,Λ−α [N] ]−1Eθ[‖θ‖2,Λ−β [p] ]2) 1 2(ζ+β−α) , (23)
in which case we have that
Eδα,β(P, p,N) . Eθ[‖θ‖22,Λ−β [p]
] ζ−α
(ζ+β−α)Eδ[‖δ‖22,Λ−α [N]
] β (ζ+β−α) . (24)
On the philosophical level, the result says that when the learning model is smoothing, that is, when the singular values of Ψ decays fast, we can select the appropriate weight matrix to compensate the smoothing effect so that the generalization error is optimized. The assumption that we have access to the exact form of the singular value decomposition of the feature matrix is only made trivialize the calculations, and is by no means essential. When the optimization algorithm is stopped before perfect matching can be achieve, the fitting from the weighted optimization scheme generated a smoother approximation (with a smaller p, according to (23), than what we would obtain with a regular least-squares minimization).
Feature matrices with fast-decaying singular values are ubiquitous in applications. In Appendix B, we provide some discussions on the applicability of this result in understanding general kernel learning (Amari & Wu, 1999; Kamnitsas et al., 2018; Jean et al., 2018; Owhadi & Yoo, 2019; Bordelon et al., 2020; Canatar et al., 2021) and learning with simplified neural networks.

heading: 6 CONCLUDING REMARKS
text: In this work, we analyzed the impact of weighted optimization on the generalization capability of feature regression with noisy data for the RFF model and generalized the result to the case of feature regression and kernel regression. For the RFF model, we show that the proposed weighting scheme (7) allows us to minimize the impact of noise in the training data while emphasizing the corresponding features according to the a priori knowledge we have on the distribution of the features. In general, emphasizing low-frequency features (i.e., searching for smoother functions) provide better generalization ability.
While what we analyze in this paper is mainly motivated by the machine learning literature, the problem of fitting models such as the RFF model (1) to the observed data is a standard inverse problem that has been extensively studied (Engl et al., 1996; Tarantola, 2005). The main focus of classical inversion theory is on the case when Ψ[N×p] is at least rank p so that there is a unique leastsquares solution to the equation Ψθ = y for any given dataset y. This corresponds to the learning problem we described above in the underparameterized regime.
In general, weighting the least-squares allows us to have an optimization algorithm that prioritize the modes that we are interested in during the iteration process. This can be seen as a preconditioning strategy from the computational optimization perspective.
It would be of great interest to derive a rigorous theory for the weighted optimization (7) for general non-convex problems (note that while the RFF model itself is nonlinear from input x to output fθ(x), the regression problem is linear). Take a general nonlinear model of the form
F (x;θ) = y,
where F is nonlinear with respect to θ. For instance, F could be a deep neural network with θ representing the parameters of the network (for instance, the weight matrices and bias vectors at different layers). We formulate the learning problem with a weighted least-squares as
θ̂ = Λ−β[p] ŵ, with ŵ = arg min w
‖Λ−α[N ] ( F (x; Λ−β[p] w)− y ) ‖22 ,
where the weight matrices are used to control the smoothness of the gradient of the optimization procedure. The linearized problem for the learning of w gives a problem of the form Ψθ̃ = ỹ with
Ψ = (Λs[N ]F (x; Λ β [p]w)) ∗(Λβ[p]) ∗(F ′)∗(x; Λβ[p]w)(Λ β [N ]) ∗, ỹ = (Λβ[p]) ∗(F ′)∗(x; Λβ[p]w)(Λ s [N ]) ∗y.
For many practical applications in computational learning and inversion, F ′ (with respect to θ) has the properties we need for Theorem 5.1 to hold. Therefore, a local theory could be obtained. The question is, can we show that the accumulation of the results through an iterative optimization procedure (e.g., a stochastic gradient descent algorithm) does not destroy the local theory so that the conclusions we have in this work would hold globally? We believe a thorough analysis along the lines of the recent work of Ma & Ying (2021) is possible with reasonable assumptions on the convergence properties of the iterative process.

heading: SUPPLEMENTARY MATERIAL
text: Supplementary material for the paper “A Generalized Weighted Optimization Method for Computational Learning and Inversion” is organized as follows.

heading: A PROOF OF MAIN RESULTS
text: We provide here the proofs for all the results that we summarized in the main part of the paper. To simplify the presentation, we introduce the following notations. For the feature matrix Ψ ∈ CN×P , we denote by
ΨT ∈ CN×p and ΨTc ∈ CN×(P−p)
the submatrices of Ψ corresponding to the first p columns and the last P − p columns respectively. In the overparameterized regime (12), we will also use the following submatrices
Ψ[N ] ∈ CN×N and Ψ[N ]\T ∈ CN×(N−p), corresponding to the first N columns of Ψ and the last N − p columns of Ψ[N ], respectively. We will use Ψ∗ to denote the Hermitian transpose of Ψ. For a P × P diagonal matrix Λ, ΛT (≡ Λ[p]) and ΛTc denote respectively the diagonal matrices of sizes p×p and (P −p)× (P −p) that contain the first p and the last P − p diagonal elements of Λ. Matrix Λ[N ] is of size N × N and used to denote the first N diagonal elements of Λ. For any column vector θ ∈ CP×1, θT and θTc denote respectively the column vectors that contain the first p elements and the last P − p elements of θ.
To study the impact of the weight matrices Λ−α[N ] and Λ −β [p] , we introduce the re-scaled version of the feature matrix Ψ, denoted by Φ, as
Φ := Λ−α[N ]ΨΛ −β [P ] . (25)
In terms of the notations above, we have that
ΦT := Λ −α [N ]ΨTΛ −β T , and ΦTc := Λ −α [N ]ΨTcΛ −β Tc . (26)
For a given column vector θ ∈ Cd and a real-valued square-matrix X ∈ Rd×d, we denote by ‖θ‖2,X := ‖Xθ‖2 = √ (Xθ)∗Xθ
the X-weighted 2-norm of θ. We will not differentiate between finite- and infinite-dimensional vectors. In the infinite-dimensional case, we simply understand the 2-norm as the usual `2-norm.
We first derive a general form for the generalization error for training with θ sampled from a distribution with given diagonal covariance matrix. It is the starting point for most of the calculations in this paper. The calculation procedure is similar to that of Lemma 1 of Xie et al. (2020). However, our result is for the case with the additional weight matrix Λ−α[N ].

heading: A.1 PROOFS OF LEMMA A.1-LEMMA A.3
text: Lemma A.1 (General form of E0α,β(P, p,N)). Let K ∈ CP×P be a diagonal matrix, and θ be drawn from a distribution such that
Eθ[θ] = 0, Eθ[θθ∗] = K . (27)
Then the generalization error for training with weighted optimization (7) can be written as E0α,β(P, p,N) = tr ( K ) + Pα,β +Qα,β , (28)
where Pα,β = tr ( Φ+T ΦTΛ −2β T Φ + T ΦTΛ β TKTΛ β T ) − 2tr ( KTΦ + T ΦT ) ,
and Qα,β = tr ( (Φ+T ) ∗Λ−2βT Φ + T ΦTcΛ β TcKTcΛ β TcΦ ∗ Tc ) , with ΦT and ΦTc given in (26).
Proof. We introduce the new variables
w = Λβ[P ]θ, and z = Λ −α [N ] y .
We can then write the solution to the weighted least-square problem (7) as
ŵT = Φ + T z, and ŵTc = 0 ,
where Φ+T is the Moore–Penrose inverse of ΦT. Using the fact that the data y contain no random noise, we write
y = ΨTΛ −β T wT + ΨTcΛ −β Tc wTc , and z = Λ −α [N ]y = ΦTwT + ΦTcwTc .
Therefore, we have, following simple linear algebra, that
‖θ̂ − θ‖22 = ‖Λ −β T (ŵT −wT)‖ 2 2 + ‖Λ −β Tc (ŵTc −wTc)‖ 2 2
= ‖Λ−βT Φ + T (ΦTwT + ΦTcwTc)−Λ −β T wT‖ 2 2 + ‖Λ −β Tc wTc‖ 2 2
= ‖Λ−βT Φ + T ΦTcwTc −Λ −β T (IT − Φ + T ΦT)wT‖ 2 2 + ‖Λ −β Tc wTc‖ 2 2 . (29)
Next, we make the following expansions:
‖Λ−βT Φ + T ΦTcwTc −Λ −β T (I− Φ + T ΦT)wT‖ 2 2 = ‖Λ −β T Φ + T ΦTcwTc‖ 2 2 + ‖Λ −β T (I− Φ + T ΦT)wT‖ 2 2 − T1,
‖Λ−βT (I− Φ + T ΦT)wT‖ 2 2 = ‖Λ −β T wT‖ 2 2 + ‖Λ −β T Φ + T ΦTwT‖ 2 2 − T2
with T1 and T2 given respectively as T1 := 2< (( Λ−βT Φ + T ΦTcwTc )∗ Λ−βT (I− Φ + T ΦT)wT ) ,
and T2 := 2< ( w∗TΛ −β T Λ −β T Φ + T ΦTwT ) .
We therefore conclude from these expansions and (29), using the linearity property of the expectation over θ, that
Eθ[‖θ̂ − θ‖22] = Eθ[‖Λ −β T wT‖ 2 2] + Eθ[‖Λ −β Tc wTc‖ 2 2]
+ Eθ[‖Λ−βT Φ + T ΦTwT‖ 2 2] + Eθ[‖Λ −β T Φ + T ΦTcwTc‖ 2 2]− Eθ[T1]− Eθ[T2] . (30)
We first observe that the first two terms in the error are simply tr ( K ) ; that is,
Eθ[‖Λ−βT wT‖ 2 2] + Eθ[‖Λ −β Tc wTc‖ 2 2] = tr
( K ) . (31)
By the assumption (27), we also have that Eθ[θTcXθ∗T] = 0 for any matrix X such that the product is well-defined. Using the relation between θ and w, we conclude that
Eθ[T1] = 0 . (32)
We then verify, using simple trace tricks, that Eθ[T2] = 2Eθ [ < ( tr ( w∗TΛ −β T Λ −β T Φ + T ΦTwT ))] = 2< ( tr ( KTΦ + T ΦT )) = 2tr ( KTΦ + T ΦT ) , (33)
where we have also used the fact that KTΛ −β T = Λ −β T KT since both matrices are diagonal, and the last step comes from the fact that Φ+T ΦT is Hermitian.
Next, we observe, using standard trace tricks and the fact that Φ+T ΦT is Hermitian, that
Eθ[‖Λ−βT Φ + T ΦTwT‖ 2 2] = Eθ[tr ( w∗TΦ + T ΦTΛ −2β T Φ + T ΦTwT ) ]
= tr ( Φ+T ΦTΛ −2β T Φ + T ΦTEθ[wTw ∗ T] ) = tr ( Φ+T ΦTΛ −2β T Φ + T ΦTΛ β TKTΛ β T ) , (34)
and
Eθ[‖Λ−βT Φ + T ΦTcwTc‖ 2 2] = Eθ[tr ( w∗TcΦ ∗ Tc(Φ + T ) ∗Λ−2βT Φ + T ΦTcwTc ) ]
= tr ( Φ∗Tc(Φ + T ) ∗Λ−2βT Φ + T ΦTcEθ[wTcw ∗ Tc ] ) = tr ( Φ∗Tc(Φ + T ) ∗Λ−2βT Φ + T ΦTcΛ β TcKTcΛ β Tc )
= tr ( (Φ+T ) ∗Λ−2βT Φ + T ΦTcΛ β TcKTcΛ β TcΦ ∗ Tc ) . (35)
The proof is complete when we put (31), (32), (33), (34) and (35) into (30).
Next we derive the general formula for the generalization error for training with noisy data. Lemma A.2 (General form of Eδα,β(P, p,N)). Let θ be given as in Lemma A.1. Then, under the assumptions (A-I)-(A-II), the generalization error for training with weighted optimization (7) from noisy data yδ can be written as
Eδα,β(P, p,N) = E0α,β(P, p,N) + Enoise(P, p,N) (36)
where E0α,β(P, p,N) is given as in (28) and Enoise(P, p,N) is given as Enoise(P, p,N) = tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Eδ[δδ ∗] ) .
The variance of the generalization error with respect to the random noise is
Varδ ( Eθ[‖θ̂ δ − θ‖22] ) = Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δδ ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]− ( Enoise(P, p,N) )2 .
(37)
Proof. We start with the following standard error decomposition
‖θ̂ δ − θ‖22 = ‖θ̂ δ − θ̂ + θ̂ − θ‖22 = ‖θ̂ − θ‖22 + ‖θ̂ δ − θ̂‖22 + 2< ( (θ̂ δ − θ̂)∗(θ̂ − θ) ) .
Using the fact that
θ̂ δ
Tc = θ̂Tc = 0 ,
the error can be simplified to
‖θ̂ δ − θ‖22 = ‖θ̂ − θ‖22 + ‖θ̂ δ
T − θ̂T‖22 + T3, (38) where
T3 = 2< ( (θ̂ δ T − θ̂T)∗(θ̂T − θT) ) .
Taking expectation with respect to θ and then the noise δ, we have that
Eδ,θ[‖θ̂ δ − θ‖22] = Eθ[‖θ̂ − θ‖22] + Eδ[‖θ̂ δ T − θ̂T‖22] + Eδ,θ[T3] . (39)
The first term on the right-hand side is simply E0α,β(P, p,N) and does not depend on δ. To evaluate the second term which does not depend on θ, we use the fact that
θ̂ δ T − θ̂T = Λ −β T Φ + T (z δ − z) = Λ−βT Φ + T Λ −α [N ](y δ − y) = Λ−βT Φ + T Λ −α [N ]δ .
This leads to
Eδ[‖θ̂ δ T − θ̂T‖22] = Eδ[‖Λ −β T Φ + T Λ −α [N ]δ‖ 2 2] = Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]
= tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Eδ[δδ ∗] ) .
The last step is to realize that Eδ,θ[T3] = 2< ( Eδ,θ [( Λ−βT Φ + T Λ −α [N ]δ )∗ (θ̂T − θT) ]) = 0
due to the assumption that δ and θ are independent.
Utilizing the formulas in (38) and (39), and the fact that Eθ[T3] = 0, we can write down the variance of the generalization error with respect to the random noise as
Varδ ( Eθ[‖θ̂ δ − θ‖22] ) = Eδ[‖θ̂ δ T − θ̂T‖42]− ( Eδ[‖θ̂ δ T − θ̂T‖22] )2
= Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δδ ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]− ( Enoise(P, p,N) )2 .
The proof is complete.
We also need the following properties on the the feature matrix Ψ of the random Fourier model. Lemma A.3. (i) For any ζ ≥ 0, we define, in the overparameterized regime (12), the matrices Π1 ∈ CN×N and Π2 ∈ CN×N respectively as
Π1 := ΨTΛ −ζ T Ψ ∗ T and Π2 = ΨTcΛ −ζ Tc Ψ ∗ Tc .
Then Π1 and Π2 admit the decomposition Π1 = Ψ[N ]ΛΠ1,ζΨ ∗ [N ] and Π2 = Ψ[N ]ΛΠ2,ζΨ ∗ [N ] respectively with ΛΠ1,ζ and ΛΠ2,ζ diagonal matrices whose m-th diagonal elements are respectively
λ (m) Π1,ζ = N−1∑ k=0 ( ν−1∑ η=0 t−ζk+Nη ) e (N) m,k, 0 ≤ m ≤ N − 1 ,
and
λ (m) Π2,ζ = N−1∑ k=0 ( µ−1∑ η=ν t−ζk+Nη ) e (N) m,k, 0 ≤ m ≤ N − 1 ,
where e(N)m,k is defined as, denoting ωN := e 2πi N ,
e (N) m,k :=
1
N N−1∑ j=0 ω (m−k)j N = { 1, k = m 0, k 6= m , 0 ≤ m, k ≤ N − 1 .
(ii) For any α ≥ 0, we define, in the underparameterized regime (13) with p < N , the matrix Ξ ∈ Cp×(N−p) as Ξ := ( Ψ∗TΛ −2α [N ] ΨT )−1 Ψ∗TΛ −2α [N ] Ψ[N ]\T .

heading: Then Ξ∗Ξ has the following representation:
text: Ξ∗Ξ = N(Ψ∗[N ]\TΣ 2α [N ]Ψ[N ]\T) −∗Ψ∗[N ]\TΛ 4αΨ[N ]\T(Ψ ∗ [N ]\TΣ 2α [N ]Ψ[N ]\T) −1 − I[N−p] .
Proof. Part (i) is simply Lemma 2 of Xie et al. (2020). It was proved by verifying that Π1 and Π2 are both circulant matrices whose eigendecompositions are standard results in linear algebra.
To prove part (ii), let us introduce the matrices P, Q, R through the relation
Ψ∗[N ]Λ −2α [N ] Ψ[N ] =
( Ψ∗TΛ −2α [N ] ΨT Ψ ∗ TΛ −2α [N ] Ψ[N ]\T
Ψ∗[N ]\TΛ −2α [N ] ΨT Ψ ∗ [N ]\TΛ −2α [N ] Ψ[N ]\T
) ≡ (
P R R∗ Q
) .
Then by standard linear algebra, we have that( Ψ∗[N ]Λ −2α [N ] Ψ[N ] )−1 = ( P−1 + P−1R(Q−R∗P−1R)−1R∗P−1 −P−1R(Q−R∗P−1R)−1 −(Q−R∗P−1R)−1R∗P−1 (Q−R∗P−1R)−1 )
since P = Ψ∗TΛ −2α [N ] ΨT is invertible.
Meanwhile, since 1√ N Ψ[N ] is unitary, we have that( Ψ∗[N ]Λ −2α [N ] Ψ[N ] )−1 = 1
N2 Ψ∗[N ]Λ 2α [N ]Ψ[N ]
= 1
N2
( Ψ∗TΛ 2α [N ]ΨT Ψ ∗ TΛ 2α [N ]Ψ[N ]\T
Ψ∗[N ]\TΛ 2α [N ]ΨT Ψ ∗ [N ]\TΛ 2α [N ]Ψ[N ]\T
) ≡ ( P̃ R̃
R̃∗ Q̃
) .
Comparing the two inverse matrices lead us to the following identity
P−1 −P−1RR̃∗ = P̃ .
This gives us that P−1R = P̃R(I[N−p] − R̃∗R)−1 .
Therefore, we have
R∗P−∗P−1R = (I[N−p] − R̃∗R)−∗R∗P̃∗P̃R(I[N−p] − R̃∗R)−1 .
Utilizing the fact that
R̃∗R = 1
N2 Ψ∗[N ]\TΛ 2α [N ]ΨTΨ ∗ TΛ −2α [N ] Ψ[N ]\T
= 1
N2 Ψ∗[N ]\TΛ 2α [N ]
( NI[N ] −Ψ[N ]\TΨ∗[N ]\T ) Λ−2α[N ] Ψ[N ]\T = I[N−p] − Q̃Q ,
we can now have the following identity
R∗P−∗P−1R = (Q̃Q)−∗R∗P̃∗P̃R(Q̃Q)−1 .
Using the formulas for P̃ and R, as well as ΨTΨ∗T = NI[N ] −Ψ[N ]\TΨ∗[N ]\T, we have
R∗P̃∗P̃R = 1
N4 Ψ∗[N ]\TΛ −2α [N ] ΨTΨ ∗ TΛ 2α [N ]ΨTΨ ∗ TΛ 2α [N ]ΨTΨ ∗ TΛ −2α [N ] Ψ[N ]\T
= 1
N3 QΨ∗[N ]\TΛ 4α [N ]Ψ[N ]\TQ−QQ̃Q̃Q .
This finally gives us
R∗P−∗P−1R = 1
N3 Q̃−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TQ̃ −1 − I[N−p] .
The proof is complete when we insert the definion of Q̃ into the result.

heading: A.2 PROOF OF THEOREM 3.1
text: We provide the proof of Theorem 3.1 here. We split the proof into the overparameterized and the underparameterized regimes.
Proof of Theorem 3.1 (Overparameterized Regime). In the overparameterized regime (12), we have that the Moore–Penrose inverse Φ+T = Φ ∗ T(ΦTΦ ∗ T) −1. Using the definitions of ΦT in (26), we can verify that Φ+T ΦTΛ −β T = Λ −β T Ψ ∗ T(ΨTΛ −2β T Ψ ∗)−1ΨTΛ −2β T ,
Λ−βT Φ + T ΦT = Λ −2β T Ψ ∗ T(ΨTΛ −2β T Ψ ∗)−1ΨTΛ −β T ,
and Φ+T ΦTΛ −2β T Φ + T ΦT = Φ + T ΦTΛ −β T Ψ∗TΨT N Λ−βT Φ + T ΦT = Λ −2β T . We therefore have tr ( Φ+T ΦTΛ −2β T Φ + T ΦTΛ β TKTΛ β T ) = tr ( KT ) . (40)
Meanwhile, a trace trick leads to tr ( KTΦ + T ΦT ) = tr ( KTΦ ∗ T(ΦTΦ ∗ T) −1ΦT ) = tr ( (ΦTΦ ∗ T) −1ΦTKTΦ ∗ T )
= tr ( (Λ−α[N ]ΨTΛ −2β T Ψ ∗ TΛ −α [N ]) −1Λ−α[N ]ΨTΛ −β T KTΛ −β T Ψ ∗ TΛ −α [N ] ) = tr ( (ΨTΛ −2β T Ψ ∗ T) −1ΨTΛ −β T KTΛ −β T Ψ ∗ T ) . (41)
Therefore, based on Lemma A.1, we have finally that Pα,β = tr ( KT ) − 2tr ( (ΨTΛ −2β T Ψ ∗ T) −1ΨTΛ −β T KTΛ −β T Ψ ∗ T ) . (42)
In a similar manner, we can check that
Qα,β = tr ( (Φ+T ) ∗Λ−2βT Φ + T ΦTcΛ β TcKTcΛ β TcΦ ∗ Tc )
= tr ( (ΦTΦ ∗ T) −∗ΦTΛ −2β T Φ ∗ T(ΦTΦ ∗ T) −1ΦTcΛ β TcKTcΛ β TcΦ ∗ Tc )
= tr ( (ΨTΛ −2β T Ψ ∗ T) −∗ΨTΛ −4β T Ψ ∗ T(ΨTΛ −2β T Ψ ∗ T) −1ΨTcKTcΨ ∗ Tc ) . (43)
We are now ready to evaluate the terms in the generalization error. First, we have that since K = cγΛ −2γ [P ] , we have, using the definition of cγ , that
tr ( K )
= cγ P−1∑ j=0 t−2γj = 1 . (44)
Second, using the results in part (i) of Lemma A.3 and the fact that 1√ N Ψ[N ] is unitary, we have
tr ( (ΨTΛ −2β T Ψ ∗ T) −1ΨTΛ −β T KTΛ −β T Ψ ∗ T ) = cγtr ( 1 N2 Ψ[N ]Λ −1 Π1,2β Ψ∗[N ]Ψ[N ]ΛΠ1,2β+2γΨ ∗ [N ] ) = cγtr ( 1 N Ψ[N ]Λ −1 Π1,2β ΛΠ1,2β+2γΨ ∗ [N ] ) = cγtr ( Λ−1Π1,2βΛΠ1,2β+2γ
) = cγ
N−1∑ m=0
∑N−1 k=0 (∑ν−1 η=0 t −2β−2γ k+Nη ) e
(N) m,k∑N−1
k=0 (∑ν−1 η=0 t −2β k+Nη ) e (N) m,k = cγ N−1∑ k=0
∑ν−1 η=0 t
−2β−2γ k+Nη∑ν−1
η=0 t −2β k+Nη
.
Therefore, Pα,β can be computed as
Pα,β = cγ p−1∑ j=0 t−2γj − 2cγ N−1∑ k=0
∑ν−1 η=0 t
−2β−2γ k+Nη∑ν−1
η=0 t −2β k+Nη
. (45)
The results in Lemma A.3 also give that
Qα,β = tr ( 1 N2 Ψ[N ]Λ −1 Π1,2β Ψ∗[N ]Ψ[N ]ΛΠ1,4βΨ ∗ [N ] 1 N2 Ψ[N ]Λ −1 Π1,2β Ψ∗[N ]Ψ[N ]ΛΠ2,2γΨ ∗ [N ] ) = tr ( Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β ΛΠ2,2γ ) ,
which means
Qα,β = N−1∑ m=0
[∑N−1 k=0 (∑ν−1 η=0 t −4β k+Nη ) e (N) m,k ][∑N−1 k=0 (∑µ−1 η=ν t −2γ k+Nη ) e (N) m,k ] [∑N−1
k=0 (∑ν−1 η=0 t −2β k+Nη ) e (N) m,k ]2 =
N−1∑ k=0
(∑ν−1 η=0 t −4β k+Nη )(∑µ−1 η=ν t −2γ k+Nη ) [∑ν−1
η=0 t −2β k+Nη ]2 . (46) We can now put (44), (45) and (46) together into the general error formula (28) to get the result (14) in Theorem 3.1. The proof is complete.
Proof of Theorem 3.1 (Underparameterized Regime). In the underparameterized regime (13), we have that the Moore–Penrose inverse Φ+T = (Φ ∗ TΦT) −1Φ∗T. This leads to
Φ+T ΦT = IT and (Φ + T ) ∗Λ−2βT Φ + T = ΦT(Φ ∗ TΦT) −∗Λ−2βT (Φ ∗ TΦT) −1Φ∗T .
Therefore, using notations in Lemma A.1, we have that Pα,β is simplified to Pα,β = −tr ( KT ) ,
while Qα,β is simplified to Qα,β = tr ( ΦT(Φ ∗ TΦT) −∗Λ−2βT (Φ ∗ TΦT) −1Φ∗TΛ −α [N ]ΨTcKTcΨ ∗ TcΛ −α [N ] ) = tr ( Λ−2α[N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] ΨTcKTcΨ ∗ Tc ) . (47)
In this regime, p ≤ N ≤ P . We therefore have Tc = ( [N ]\T ) ∪ [N ]c ([N ]c := [P ]\[N ]). Using the fact that K is diagonal, we obtain
ΨTcKTcΨ ∗ Tc = Ψ[N ]\TK[N ]\TΨ ∗ [N ]\T + Ψ[N ]cK[N ]cΨ ∗ [N ]c
Following the result of Lemma A.3, the second part of decomposition is simply the matrix Π2 with ν = 1. To avoid confusion, we denote it by Π3 and use the result of Lemma A.3 to get
Ψ[N ]cK[N ]cΨ ∗ [N ]c = Ψ[N ]ΛΠ3,2γΨ ∗ [N ]
where the diagonal elements of ΛΠ3,2γ are
λ (m) Π3,2γ = N−1∑ k=0 ( µ−1∑ η=1 t−2γk+Nη )( 1 N N−1∑ j=0 ω (m−k)j N ) , 0 ≤ m ≤ N − 1 . (48)
Next we perform the decomposition
Ψ[N ]ΛΠ3,2γΨ ∗ [N ] = ΨTΛΠ3,2γTΨ ∗ T + Ψ[N ]\TΛΠ3,2γ [N ]\TΨ ∗ [N ]\T ,
and define the diagonal matrices
K̂1 := ΛΠ3,2γT and K̂2 := K[N ]\T + ΛΠ3,2γ [N ]\T .
We can then have ΨTcKTcΨ ∗ Tc = ΨTK̂1Ψ ∗ T + Ψ[N ]\TK̂2Ψ ∗ [N ]\T .
Plugging this into the expression for Qα,β , we have Qα,β = tr ( Λ−2α[N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] ΨTK̂1Ψ ∗ T )
+ tr ( Λ−2α[N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\TK̂2Ψ ∗ [N ]\T ) . (49)
The first term simplifies to tr ( K̂1 ) . The second term vanishes in the case of α = 0 and in the case when the problem is formally determined; that is, when p = N , but is in general nonzero when p < N . Using the result in part (ii) of Lemma A.3, we have that
tr ( Λ−2α[N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\TK̂2Ψ ∗ [N ]\T ) = tr ( Ψ∗[N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\TK̂2
) = Ntr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1K̂2 ) − tr ( K̂2 ) , X := Ψ∗[N ]\TΛ 2α [N ]Ψ[N ]\T .
Using this result and the singular value decomposition
Σα[N ]Ψ[N ]\T = Udiag([Σ00, · · · ,Σ(N−p−1)(N−p−1)])V∗
introduced in Theorem 3.1, we have Qα,β = tr ( K̂1 − K̂2 ) +Ntr ( Σ−1U∗Λ2α[N ]UΣ −1VK̂2V ∗)
= tr ( ΛΠ3,2γT −K[N ]\T −ΛΠ3,2γ [N ]\T ) +N N−p−1∑ i=0 N−p−1∑ j=0 ẽ (N) ij ê (N) ji ΣiiΣjj .
where
ẽ (N) ij = N−1∑ k=0 t2αk UkiUkj , ê (N) ij = N−p−1∑ k′=0 (cγt −2γ p+k′ + λ (p+k′) Π3,2γ )V ik′Vjk′ , 0 ≤ i, j ≤ N − p− 1 .
The proof is complete when we insert the expressions of Pα,β and Qα,β back to the general formula (28) and use the form of λ(m)Π3,2γ in (48).
When the problem is formally determined, i.e., in the case of N = p (equivalent to ν = 1), Ψ∗TΛ −2α [N ] ΨT = ΨTΛ −2α [N ] Ψ ∗ T. This allows us to find that
Qα,β = tr ( 1 N2 Ψ[N ]Λ −1 Π1,2α Ψ∗[N ]Ψ[N ]ΛΠ1,4αΨ ∗ [N ]Ψ[N ]Λ −1 Π1,2α Ψ∗[N ] 1 N2 Ψ[N ]ΛΠ2,2γΨ ∗ [N ] ) = tr ( Λ−2Π1,2αΛΠ1,4αΛΠ2,2γ ) =
N−1∑ k=0
(∑ν−1 η=0 t −4α k+Nη )(∑µ−1 η=ν t −2γ k+Nη ) [∑ν−1
η=0 t −2α k+Nη ]2 , which degenerates to its form (46) in the overparameterized regime with ν = 1. The result is then independent of α since the terms that involve α cancel each other. Similarly, the simplification of Qα,β in (46) with ν = 1 (N = p) also leads to the fact that β disappears from the formula.

heading: A.3 PROOFS OF LEMMA 4.1 AND THEOREM 4.2
text: Proof of Lemma 4.1. By Lemma A.2, the main task is to estimate the size of the generalization error caused by the random noise in the two regimes. The error is, according to (36),
Enoise(P, p,N) = tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Eθ[δδ ∗] ) .
In the overparameterized regime, we have that
Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] = Λ −α [N ](ΦTΦ ∗ T) −∗ΦTΛ −2β T Φ ∗ T(ΦTΦ ∗ T) −1Λ−α[N ]
= (ΨTΛ −2β T Ψ ∗ T) −∗ΨTΛ −4β T Ψ ∗ T(ΨTΛ −2β T Ψ ∗ T) −1 .
Therefore we have, using the results in part (i) of Lemma A.3,
Enoise(P, p,N) = σ2tr ( 1 N2 Ψ[N ]Λ −1 Π1,2β Ψ∗[N ]Ψ[N ]ΛΠ1,4βΨ ∗ [N ] 1 N2 Ψ[N ]Λ −1 Π1,2β Ψ∗[N ] )
= σ2 N tr ( Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β ) = σ2 N N−1∑ k=0
∑ν−1 η=0 t
−4β k+Nη[∑ν−1
η=0 t −2β k+Nη ]2 . To get the variance of the generalization error with respect to noise, we first compute
Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δδ ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]
= σ4tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) +2σ4tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Λ −α [N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]
) = σ4
N4 tr ( Ψ[N ]Λ −1 Π1,2β ΛΠ1,4βΛ −1 Π1,2β Ψ∗[N ] ) tr ( Ψ[N ]Λ −1 Π1,2β ΛΠ1,4βΛ −1 Π1,2β Ψ[N ] ) + 2σ4
N4 tr ( Ψ[N ]Λ −1 Π1,2β ΛΠ1,4βΛ −1 Π1,2β Ψ∗[N ]Ψ[N ]Λ −1 Π1,2β ΛΠ1,4βΛ −1 Π1,2β Ψ∗[N ] )
= σ4 N2 tr ( Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β ) tr ( Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β ) + 2σ4
N2 tr ( Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β Λ−1Π1,2βΛΠ1,4βΛ −1 Π1,2β ) = σ4
N2 [N−1∑ k=0 ∑ν−1 η=0 t −4β k+Nη[∑ν−1
η=0 t −2β k+Nη
]2 ]2 + 2σ4N2 N−1∑ k=0
[∑ν−1 η=0 t −4β k+Nη ]2 [∑ν−1
η=0 t −2β k+Nη
]4 .
This, together with the general form of the variance in (37), gives the result in (18).
In the underparameterization regime, we have that
Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] = Λ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] .
Therefore we have, using the fact that ΨTΨ∗T + Ψ[N ]\TΨ ∗ [N ]\T = NI[N ], that
Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]
= 1
N2
( ΨTΨ ∗ T + Ψ[N ]\TΨ ∗ [N ]\T ) Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ( ΨTΨ ∗ T + Ψ[N ]\TΨ ∗ [N ]\T ) = 1
N2 ΨTΨ
∗ T +
1
N2 ΨT(Ψ
∗ TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\TΨ ∗ [N ]\T
+ 1
N2 Ψ[N ]\TΨ
∗ [N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗Ψ∗T
+ 1
N2 Ψ[N ]\TΨ
∗ [N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\TΨ ∗ [N ]\T . (50)
This leads to, using the fact that Ψ∗TΨT = NI[p] and properties of traces,
Enoise(P, p,N) = σ2tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) = p
N σ2 +
σ2 N tr ( Ψ∗[N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\T ) Using the result in part (ii) of Lemma A.3, we have that the second term simplifies to
tr ( Ψ∗[N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\T ) = Ntr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX
−1)− (N − p), X := Ψ∗[N ]\TΛ2α[N ]Ψ[N ]\T . Using this result and the singular value decomposition
Σα[N ]Ψ[N ]\T = Udiag([Σ00, · · · ,Σ(N−p−1)(N−p−1)])V∗
introduced in Theorem 3.1, we have
Enoise(P, p,N) = p N σ2−N − p N σ2+σ2tr
( Σ−1U∗Λ2α[N ]UΣ −1) = σ2(2p−N N + N−p−1∑ j=0 ẽ (N) jj Σ2jj ) .
Following the same procedure as in the overparameterization regime, the variance of the generalization error with respect to noise in this case is
Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δδ ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]
= σ4tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) + 2σ4tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Λ −α [N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] )
The first term on the right hand side is simply ( Enoise(P, p,N) )2 . To evaluate the second term, we use the formula (50). It is straightforward to obtain, after some algebra, that
tr ( Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]Λ −α [N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ] ) = 1
N3 tr ( ΨTΨ ∗ T ) + 2 N2 tr ( Ψ∗[N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\T ) + 1
N2 tr ([ Ψ∗[N ]\TΛ −2α [N ] ΨT(Ψ ∗ TΛ −2α [N ] ΨT) −∗(Ψ∗TΛ −2α [N ] ΨT) −1Ψ∗TΛ −2α [N ] Ψ[N ]\T ]2) = p
N2 +
2
N2
( Ntr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1)− (N − p)) + 1
N2
( N2tr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1) −2Ntr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX
−1)+ tr(I[N−p])) =
2p−N N2
+ tr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1) =
2p−N N2
+ tr ( X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1X−∗Ψ∗[N ]\TΛ 4α [N ]Ψ[N ]\TX −1) This leads to
Eδ[tr ( δ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δδ ∗Λ−α[N ](Φ + T ) ∗Λ−2βT Φ + T Λ −α [N ]δ ) ]
= ( Enoise(P, p,N) )2 + 2σ4 (2p−N N2 + tr ( Σ−2U∗TΛ 2α [N ]UTΣ −2U∗TΛ 2α [N ]UT )) = ( Enoise(P, p,N) )2 + 2σ4 (2p−N N2 + N−p−1∑ i=0 N−p−1∑ j=0 ẽ (N) ij ẽ (N) ji Σ2iiΣ 2 jj ) .
Inserting this into the general form of the variance, i.e. (37), will lead to (19) for the underparameterized regime. The proof is complete.
Proof of Theorem 4.2. By standard decomposition, we have that
‖θ̂ δ − θ‖22 = ‖θ̂ δ − θ̂ + θ̂ − θ‖22
= ∥∥∥∥( Ψ+TO(P−p)×N ) (yδ − y) + (I[P ] − ( Ψ+T Ψ O(P−p)×P ) )θ ∥∥∥∥2 2
≤ 2‖Ψ+T δ‖ 2 2 + 2 ∥∥∥∥(I[P ] − ( Ψ+T ΨO(P−p)×P ) )θ ∥∥∥∥2 2 . (51)
From the calculations in the previous sections, it is straightforward to verify that, under the scaling of the parameters assumed in the theorem, we have
‖Ψ+T ‖2,Λ−α [N] 7→2 ∼ p − 12−α and ∥∥∥∥I[P ] − ( Ψ+T ΨO(P−p)×P )∥∥∥∥ 2,Λ−β [P ] 7→2 ∼ p−β(1− p−2γ+1) ,
where ‖X‖2,Λ−α [N] 7→2 denote the norm of X as an operator from the Λ −α [N ]-weighted C N to Cp.
This allows us to conclude that, after taking expectation with respect to θ and then δ,
Eδ,θ[‖θ̂ δ − θ‖22] ≤ 2‖Ψ+‖22,Λ−α [N] 7→2Eδ[‖Λ −α [N ](y δ − y)‖22]
+2‖(I[P ]−Ψ+Ψ)‖22,Λ−β [P ] 7→2Eθ[‖θ‖ 2 2,Λ−β [P ] ] . p−2α−1Eδ[‖δ‖22,Λ−α [N] ]+p−2β(1−p−2γ+1)Eθ[‖θ‖22,Λ−β [P ] ].
When α > 0, the error decrease monotonically with respect to p. When α < −1/2, the first term increase with p while the second term descreases with p. To minimize the right-hand side, we take
p ∼ (Eδ[‖δ‖22,Λ−α [N] ]−1Eθ[‖θ‖22,Λ−β [P ] ]) 1 2(β−α) . (52)
This leads to Eθ,δ[‖θ̂ δ − θ‖22] . Eθ[‖θ‖22,Λ−β
[P ]
] −2α+1 2(β−α)Eδ[‖δ‖22,Λ−α
[N]
] 2β+1 2(β−α) . (53)
The proof is now complete.

heading: A.4 PROOF OF THEOREM 5.1
text: Proof of Theorem 5.1. Without loss of generality, we assume that Ψ is diagonal with diagonal elements Ψkk ∼ k−ζ . Let Ψ+ be the pseudoinverse of Ψ that consists of the first p features such that
(Ψ+)qq ∼ { qζ , q ≤ p 0, q > p
where ζ is the exponent in the SVD of Ψ assumed in Theorem 5.1.
We can then check that the operators I−Ψ+Ψ : `2 Λ−β
[P ]
7→ `2 and Ψ+ : `2 Λ−α
[N] 7→ `2 have the following norms respectively
‖Ψ+‖2,Λ−α [N] 7→2 ∼ p ζ−α and ‖(I−Ψ+Ψ)‖2,Λ−β [P ] 7→2 ∼ p −β .
By the error decomposition (51), we then conclude, after taking expectation with respect to θ and then δ, that
Eδ,θ[‖θ̂ δ
c − θ‖22] . ‖Ψ+‖22,Λ−α [N] 7→2Eδ[‖Λ −α [N ]δ‖ 2 2]
+ ‖(I−Ψ+Ψ)‖2 2,Λ−β [P ] 7→2Eθ[‖θ‖ 2 2,Λ−β [P ] ] . p2(ζ−α)Eδ[‖δ‖22,Λ−α [N] ] + p−2βEθ[‖θ‖22,Λ−β [P ] ].
We can now select p ∼ (Eδ[‖δ‖22,Λ−α
[N]
]−1Eθ[‖θ‖2,Λ−β [P ] ]2) 1 2(ζ+β−α) (54)
to minimize the error. This leads to
Eδ,θ[‖θ̂ δ − θ‖22] . Eθ[‖θ‖22,Λ−β
[P ]
] ζ−α
(ζ−α+β)Eδ[‖δ‖22,Λ−α [N]
] β (γ+β−α) . (55)
The proof is now complete.

heading: B ON THE APPLICABILITY OF THEOREM 5.1
text: In order for the result of Theorem 5.1 to hold, we need the model to be learned to have the smoothing property. In the case of feature regression, this requires that the feature matrices (or sampling matrices) corresponding to those kernels have singular values decaying algebraically. This turns out to be true for many kernel regression models in practical applications. In the case of solving inverse problems, this is extremely common as most inverse problems based on physical models have smoothing operators; see, for instance, Isakov (2006); Kirsch (2011) for a more in-depth discussion on this issue.
General kernel regression. Let H be a reproducing kernel Hilbert space (RKHS) over X , and K : H × H → R be the corresponding (symmetric and positive semidefinite) reproducing kernel. We are interested in learning a function f∗ from given data {xj , yδj}Nj=1. The learning process can be formulated as
min f∈H N∑ j=1 ( f(xj)− yδj )2 + β‖f‖2H . (56)
Let {ϕk}k≥0 be the eigenfunctions of the kernel K such that∫ K(x, x̃)ϕk(x̃)µ(x̃)dx̃ = λkϕk(x) , (57)
where µ is the probability measure that generates the data. Following Mercer’s Theorem (Rasmussen & Williams, 2006), we have thatK admits a representation in terms of its kernel eigenfunctions; that is, K(x, x̃) = ∑ k≥0 λkϕk(x)ϕk(x̃) . Moreover, the solution to the learning problem as well as the target function can be approximated respectively as
f∗θ (x) = P−1∑ k=0 θkϕk(x), and fθ̂p(x) = p−1∑ k=0 θ̂kϕk(x) , (58)
where to be consistent with the setup in the previous sections, we have used P and p to denote the numbers of modes in the target function and the learning solution respectively. We can now define Ψ to be the feature matrix (with components (Ψ)kj = ϕk(xj)) so that the kernel regression problem can be recast as the optimization problem
θ̂ δ
p = arg min θ
‖Ψθp − yδ‖22 + β‖θp‖22 . (59)
For a given training data set consisting of N data points, the generalization error for this learning problem can then be written as in the same form as (6); that is,
Eδα,β(P, p,N) = Eθ,δ [ ‖f∗θ (x)− fθ̂p(x)‖ 2 L2(X) ] = Eθ,δ [ ‖θ̂p − θ‖22 ] , (60)
using the generalized Parseval’s identity (also Mercer’s Theorem) in the corresponding reproducing kernel Hilbert space.
Popular kernels (Rasmussen & Williams, 2006, Chapter 4) in applications include the polynomial kernel KPolynomial(x, x̃) = (α〈x, x̃〉+ 1)d, (61) where d is the dimension of the underlying space, and the Gaussian RBF (Radial Basis Function) kernel
KGaussian(x, x̃) = exp
( −‖x− x̃‖ 2
2σ2
) . (62)
For appropriate datasets, such as those normalized ones that live on the unit sphere Sd−1, there is theoretical as well as numerical evidence to show that the feature matrix Ψ has eigenvalues decay fast (for instance, algebraically). This means that for such problems, we also have that the lowfrequency modes dominate the high-frequency modes in the target function. This means that the weighted optimization framework we analyzed in the previous sections should also apply here. We refer interested readers to Rasmussen & Williams (2006) and references therein for more technical details and summarize the main theoretical results in Appendix B.
Neural tangent kernel. It turns out that a similar technique can be used to understand some aspects of learning with neural networks. It is particularly related to the frequency bias of neural networks that has been extensively studied (Ronen et al., 2019; Wang et al., 2020). It is also closely related to the regularization properties of neural networks (Martin & Mahoney, 2018). To make the connection, we consider the training of a simple two-layer neural network following the work of Yang & Salman (2019) and Ronen et al. (2019). We refer interested readers to Daniely et al. (2016) where kernel formulation of the initialization of deep neural nets was first introduced. In their setting, the neural network is a concentration of a computation skeleton, and the initialization of the neural network is done by sampling a Gaussian random variable.
We denote by f a two-layer neural network that takes input vector x ∈ Rd to output a scalar value. We assume that the hidden layer has J neurons. We can then write the network, with activation function σ, as
f(x; Θ,α) = 1√ M M∑ m=1 αmσ(θ T mx), (63)
where Θ = [θ1, · · · ,θM ] ∈ Rd×M and α = [α1, · · · , αM ]T ∈ RJ are respectively the weights of the hidden layer and the output layer of the network. We omit bias in the model only for simplicity.
In the analysis of Ronen et al. (2019); Yang & Salman (2019), it is assumed that the weight of the output layer α is known and one is therefore only interested in fitting the data to get Θ. This training
process is done by minimizing the L2 loss over the data set {xj , yj}Nj=1:
Φ(Θ) = 1
2 N∑ j=1 ( yj − f(xj ; Θ,α) )2 .
When the activation function σ is taken as the ReLU function σ(x) = max(x, 0), we can define the matrix, which depends on Θ,
Ψ(Θ) = 1√ M
 α1χ11x1 α2χ12x1 · · · αMχ1Mx1 α1χ21x2 α2χ22x2 · · · αMχ2Mx2
... ...
. . . ...
α1χN1xN α2χN2xN · · · αMχNMxN  , where χjm = 1 if θTmxj ≥ 0 and χjm = 0 if θ T mxj < 0. The least-squares training loss can then be written as 12‖Ψ(Θ)Θ− y‖ 2 2. A linearization of the least-squares problem around Θ0 can then be formulated as ∆Θ = arg min
∆Θ ‖Ψ(Θ0)∆Θ−∆y‖22, (64)
where ∆y := y −Ψ(Θ0)Θ0 is the perturbed data. Under the assumption that the input data are normalized such that ‖x‖ = 1 and the weight αk ∼ U(−1, 1) (1 ≤ k ≤ M ), it was shown in Ronen et al. (2019) that the singular values of the matrix Ψ(Θ0) decays algebraically. In fact, starting from initialization θm ∼ N (0, κ2I), this result holds during the whole training process under some mild assumptions.
We summarize the main result on the Gaussian kernel and the linear kernel in the following theorem. Theorem B.1 (Theorem 2 and Theorem 3 of Minh et al. (2006)). Let X = Sd−1, d ∈ N and d ≥ 2. Let µ be the uniform probability distribution on Sd−1. Then eigenvalues and eigenfunctions for the Gaussian kernel are respectively
λk = e −2/σ2σd−2Ik+n/2−1
( 2
σ2
) Γ (n
2
) (65)
for all k ≥ 0, where I denotes the odified Bessel function of the first kind. Each λk occurs with multiplicity N(n, k) with the corresponding eigenfunctions being spherical harmonics of order k on Sd−1. The λk’s are decreasing if σ ≥ √ 2/d for the Gaussian kernel (Rasmussen & Williams, 2006). Furthermore, λk forms a descreasing sequence and( 2e
σ2
)k A1
(2k + n− 2)k+n−12 < λk
( 2e
σ2
)k A2
(2k + n− 2)k+n−12 The nonzero eigenvalues are
λk = 2 d+n−2 d! (d− k)! Γ(d+ n−12 )Γ( n 2 )√ πΓ(d+ k + n− 1) (66)
for the polynommial kernel with 0 ≤ k ≤ d. Each λk occurs with multiplicity N(n, k), with the correpsonding eigenfunctions being spherical harmonics of order k on Sn−1. Furthermore, λk forms a descreasing sequence and
B1
(k + d+ n− 2)2d+n− 32 < λk <
B2
(k + d+ n− 2)d+n− 32 .
The results have been proved in different settings. When the samples are not on Sd−1 but in the cube [−1, 1]d or the unit ball Bd, or the underlying distrubion of the data is not uniform, there are similar results. We refer interested readers to Minh et al. (2006); Rasmussen & Williams (2006) and references therein for more details. In Figure 2, we plot the singular values of the sampling matrix for the polynomial kernel and the Gaussian RBF kernel for the MINST data set.
For the theorem to work for learning with neural networks as we discussed in Section 5, it has been shown that the neural tangent kernel also satisfies the required property in different settings. The main argument is that neural tangent kernel is equivalent to kernel regression with the Laplace kernel as proved in Chen & Xu (2020). We cite the following result for the two-layer neural network model.
Theorem B.2 (Proposition 5 of Bietti & Mairal (2019)). For any x, x̃ ∈ Sd−1, the eigenvalues of the neural tangent kernel K are non-negative, satisfying µ0, µ1 > 0, µk = 0 if k = 2j + 1 with j ≥ 1, and otherwise µk ∼ C(d)k−d as k → ∞, with C(d) a constant depending only on d. The eigenfunction corresponds µk is the spherical harmonic polynomials of degree k.
The proof of this result can be found in Bietti & Mairal (2019). The similarity between neural tangent kernel and the Laplace kernel was first documented in Geifman et al. (2020), and a rigorous theory was developed in Chen & Xu (2020).

Review summaries:
This paper studies the weighted least squares, random features model under noise one-dimensional data setting in under-/over-parameterized regime. The derived error bounds demonstrate the impact of noise on the generalization error. Besides, the extension to kernel regression shows that, the selected weighted matrix is helpful to generalization when the RKHS is small (i.e., singular values of \Psi decay fast).

 The paper follows and extends the work by Belkin, Hsu, and Xu (2020) and  Liang and Rakhlin (2018) which studies the bias-variance trade-off of the regression/interpolation problem in the under/over-parametrization regions.  In particular, this paper follows the random Fourier model setting in Xie et al. (2020) and analyzed a generalized weighted least-square optimization method that allows the weighting in both the parametrization and data space. The authors derived the generalization error of such weighted least-square framework for the over parametrized and under parameterized regimes and compare them in these two cases.  The general conclusion is that emphasizing low-frequency features provide better generalization ability.  The paper also studies both noise-free and noise cases.