Can we globally optimize cross-validation loss? Quasiconvexity in ridge regression
Year: 2022

Abstract:
Models like LASSO and ridge regression are extensively used in practice due to their interpretability, ease of use, and strong theoretical guarantees. Crossvalidation (CV) is widely used for hyperparameter tuning in these models, but do practical optimization methods minimize the true out-of-sample loss? A recent line of research promises to show that the optimum of the CV loss matches the optimum of the out-of-sample loss (possibly after simple corrections). It remains to show how tractable it is to minimize the CV loss. In the present paper, we show that, in the case of ridge regression, the CV loss may fail to be quasiconvex and thus may have multiple local optima. We can guarantee that the CV loss is quasiconvex in at least one case: when the spectrum of the covariate matrix is nearly flat and the noise in the observed responses is not too high. More generally, we show that quasiconvexity status is independent of many properties of the observed data (response norm, covariate-matrix right singular vectors, and singular-value scaling) and has a complex dependence on the few that remain. We empirically confirm our theory using simulated experiments.

heading: 1 Introduction
text: Linear models, including LASSO and ridge regression, are widely used for data analysis across diverse applied disciplines. Linear models are often preferred since they are straightforward to apply in various senses. In particular, (1) their parameters are readily interpretable. (2) They have strong theoretical guarantees on quality. And (3) standard optimization tools are often assumed to find useful parameter and hyperparameter values. Despite their seeming simplicity, though, mysteries remain about the quality of inference in linear models. Consider cross-validation (CV) [Stone, 1974, Allen, 1974], the de facto standard for hyperparameter selection across machine learning methods [Musgrave et al., 2020]. CV is an easy-to-evaluate proxy for the true out-of-sample loss. Is it a good proxy? [Homrighausen and McDonald, 2014, 2013, Chetverikov et al., 2020, Hastie et al., 2020, Patil et al., 2021] give conditions under which the global minimum of the CV loss (possibly with some mild corrections) matches the optimum of the out-of-sample loss in LASSO and ridge regression. To complete the picture, we must understand whether standard methods for minimizing the CV loss find a global minimum.
It would be easy to find a unique minimum of the CV loss if the CV loss were convex. Alas (though perhaps unsurprisingly), we show below that in essentially every case of interest the CV loss is not convex. Indeed, the usual introductory cartoon of CV loss (left panel of Fig. 1; see also Fig. 5.9 of Hastie et al. [2017] or Fig. 1 of Rad and Maleki [2020]) is not convex. But the cartoon CV loss still exhibits a single global minimum and would be easy to globally minimize with popular approaches like gradient-based methods [Do et al., 2007, Maclaurin et al., 2015, Pedregosa, 2016, Lorraine et al., 2020] or grid search [Bergstra and Bengio, 2012, Pedregosa et al., 2011, Hsu et al., 2003]. Indeed, a
∗Alternate email: wtstephe@gmail.com
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
more plausible possibility (which holds for the typical cartoon CV loss) is that the CV loss might be quasiconvex; i.e. its level sets are convex. The benefit of quasiconvexity is that, in one dimension, any local optimum is a global optimum.
Unfortunately, this cartoon need not hold in general, even in simple models like ℓ2-regularized linear regression (i.e. ridge regression). Consider minimizing leave-one-out CV (LOOCV) loss as a function of the ridge regularization parameter; we denote this loss by L. Wilson et al. [2020, Fig. 1] detail a simulated example in which L can be non-quasiconvex. We first demonstrate that L can be non-quasiconvex in real-data examples; see the middle and right panel of Fig. 1, which we describe in detail in Section 3. We next characterize which aspects of the covariate matrix and observed responses affect quasiconvexity. We prove that the norm of the responses, the scale of the singular values of the covariate matrix, and the right singular vectors of the covariate matrix all have no effect on the quasiconvexity of L. While this result places substantial constraints on what drives the quasiconvexity of L, we show that the quasiconvexity of L is unfortunately still a complex function of the remaining quantities. Our third contribution is to prove conditions under which L is guaranteed to be quasiconvex. In particular, we show that if (1) the covariate matrix has a singular value spectrum sufficiently close to uniform, (2) the least-squares estimator fits the training data sufficiently well, and (3) the left singular vectors of the covariate matrix are sufficiently regular, then L is guaranteed to be quasiconvex. While the conditions of our theory are deterministic, we show that they have natural probabilistic interpretations; as a corollary to our theory, we demonstrate that many of our conditions are satisfied either empirically or theoretically by well-specified linear regression problems with i.i.d. sub-Gaussian covariates and moderate signal-to-noise ratios. Through empirical studies, we validate the conclusions of our theory and the necessity of our assumptions.

heading: 2 Setup and notation
text: For n ∈ {1, . . . , N}, we observe covariates xn ∈ RD and responses yn ∈ R. We are interested in learning a linear model between the covariates and responses, ⟨xn, θ⟩ ≈ yn, for some parameter θ ∈ RD. In ridge regression, i.e. ℓ2-regularized linear regression, we take some λ > 0 and estimate:
θ̂(λ) := argmin θ∈RD N∑ n=1 (⟨xn, θ⟩ − yn)2 + λ 2 ∥θ∥22 . (1)
The regularization parameter λ is typically chosen by minimizing the cross-validation (CV) loss. Here we study the leave-one-out CV (LOOCV) loss:
L(λ) := N∑
n=1
( ⟨xn, θ̂\n(λ)⟩ − yn )2 , (2)
where θ̂\n(λ) is the solution to Eq. (1) with the nth datapoint left out.
Let the covariate matrix X ∈ RN×D be the matrix with rows xn, and let the vector Y ∈ RN be the vector with entries yn. We consider the low to modest-dimensional case where D < N and assume the covariate matrix X is full-rank. We further assume X and Y have undergone standard data pre-processing, as described next.
Condition 1. Y is zero-mean, and X has zero-mean, unit variance columns. Equivalently, where 1 ∈ RN is the vector of all ones, 1TY = 0 and XT1 = 0 ∈ RD and for all d = 1, . . . , D,∑N
n=1 x 2 nd = N .
Preprocessing X and Y to satisfy Condition 1 represents standard best practice for ridge regression. First, using an unregularized bias parameter in Eq. (1) and setting Y to be zero-mean are equivalent; we choose to make Y zero-mean, as it simplifies our analysis below. The conditions on the covariate matrix X are important to ensure the use of ℓ2-regularization is sensible. In particular, Eq. (1) penalizes all coordinates of θ equally. If e.g. some columns of X are measured in different scales or are centered differently, this uniform penalty will be inappropriate.

heading: 3 LOOCV loss is typically not convex and need not be quasiconvex
text: If the LOOCV loss L were convex or quasiconvex in λ, then any local minimum of L would be a global minimum, and we could trust gradient-based optimization methods or grid search methods to return a value near a global minimum. We next see that unfortunately L is typically not convex and is often not even quasiconvex. First we show that, in essentially all cases of interest, L is not convex. Proposition 1. If λ = ∞ is not a minimum of L, then L is not a convex function.
Proof. For the sake of contradiction, assume L is convex and λ = ∞ is not a minimum of L. This implies that there is some maximal λ∗ < ∞ such that L′(λ∗) = 0. Let δ := L′(λ∗ + 1). By convexity, L′′ ≥ 0, so we know that δ > 0 and that for λ ≥ λ∗ + 1, we have L′(λ) ≥ δ. Thus for λ ≥ λ∗ + 1, we have L(λ) ≥ δ(λ − λ∗ − 1). So limλ→∞ L(λ) = ∞. However, inspection of L shows limλ→∞ L(λ) = ∑N n=1 y 2 n < ∞, which is a contradiction.
We say that the result covers essentially all cases of interest: if L continues to decrease as λ → ∞, then there is so little signal in the data that the zero model θ = 0 ∈ RD is the optimal predictor according to LOOCV.
Although L is generally not convex, L(λ) might still be easy to optimize if it satisfies an appropriate generalized notion of convexity. To that end, we recall the definition of quasiconvexity.
Definition 1. A function f : Rp → R is quasiconvex if its level sets are convex.
In one dimension (i.e. p = 1 in Definition 1), quasiconvexity guarantees that any local optimum is a global optimum, just as convexity does. This property is often the key consideration in practical optimization algorithms. Moreover, it is not unreasonable to hope that the CV loss is quasiconvex: typical illustrations of the CV loss are not convex but are quasiconvex; see e.g. Hastie et al. [2015,
Fig. 5.9], Rad and Maleki [2020, Fig. 1], or the left panel of Fig. 1. Illustrations of the out-of-sample loss are also typically quasiconvex; see e.g. Fig. 3.6 of Bishop [2006].
Unfortunately, we next demonstrate that the CV loss derived from real data analysis problems can be non-quasiconvex. Our first dataset contains N = 2,938 observations of life expectancy, along with D = 20 covariates such as country of origin or alcohol use; see ?? for a full description. In this case, after pre-processing according to Condition 1, L for the full dataset is quasiconvex. But now consider some additional standard data pre-processing. Practitioners often perform principal component regression (PCR) with the aim of reducing noise in the estimated θ. That is, they take the singular value decomposition of X = USV and produce an N × R dimensional covariate matrix X ′ by retaining just the top R singular values of X: X ′ = U·,:RS:R. With this pre-processing, the resulting LOOCV curve L is non-quasiconvex for many values of R; in the center panel of Fig. 1 we show one example for R = 15.
Our second dataset consists of recorded wine quality of N = 1,599 red wines. The goal is to predict wine quality from D = 11 observed covariates relating to the chemical properties of each wine; see ?? for a full description. We find that subsets of this dataset often exhibit non-quasiconvex L. In the right panel of Fig. 1, we show L for a subset of size N = 50. We see that this plot contains at least two local optima, with substantially different values of λ and substantially different values of the loss. A gradient-based algorithm initialized sufficiently far to the left would not find the global optimum, and grid search without sufficiently large values would not find the global optimum.
Now we know that L can be non-quasiconvex for real data. Given the difficulty of optimizing a function with several local minima, we next seek to understand when L is quasiconvex or not.

heading: 4 What does the quasiconvexity of L depend on?
text: We have seen that L can be quasiconvex or non-quasiconvex, depending on the data at hand. If we could determine the quasiconvexity of L from the data, we might better understand how to efficiently tune hyperparameters from the CV loss. In what follows, we start by showing that the quasiconvexity of L is, in fact, independent of many aspects of the data (Proposition 2). We will see, however, that the dependence of quasiconvexity on the remaining aspects (though they are few) is complex.
A linear regression problem has a number of moving parts. The response vector Y may be an arbitrary vector in RN , and the covariate matrix X can be written in terms of its singular values and left and right singular vectors. More precisely, let X = U diag(S)V T be the singular value decomposition of the covariate matrix X , where U ∈ RN×D is an N ×D matrix with orthonormal columns, S ∈ RD is a vector with positive entries, and V ∈ RD×D is an orthonormal matrix. Note we use the “compact” singular value decomposition, where U is a N ×D matrix, rather than a full N ×N matrix. With this notation in hand, we can identify aspects of the problem that do not contribute to quasiconvexity in the following result, which is proved in ??. Proposition 2. The quasiconvexity of L is independent of
1. the matrix of right singular vectors, V ,
2. the norm of the responses, ∥Y ∥2, and
3. the scaling of the singular values (i.e. changing S into S/c for c ∈ R>0),
in the sense that altering any of these quantities does not change whether or not L is quasiconvex. Remark 1. Assume Condition 1 holds. Then by Proposition 2, without loss of generality we may (and do) assume that V = ID and that Y is a vector on the unit (N − 1)-sphere.
Recall X has zero-mean columns by pre-processing the data (Condition 1). By Proposition 2, we assume without loss of generality V = ID. Thus, the columns of X have zero mean when UT1 = 0, where 1 ∈ RN is the vector of all ones. Also, while Remark 1 notes that we can consider Y to be on the unit (N − 1)-sphere, note that the condition 1TY = 0 from Condition 1 allows us to further constrain Y to be parameterized by a vector on the unit (N − 2)-sphere. Hence the quasiconvexity of L depends on three quantities: (1) the matrix of left singular vectors, U , an orthonormal matrix with 1 ∈ RN in its left null-space, (2) the (normalized) vector Y which sits on the unit (N − 2)-sphere, and (3) the (normalized) singular values.
Now that we know the quasiconvexity of L depends on only three quantities, we might hope that quasiconvexity would be a simple function of the three. To investigate this dependence, we consider the case of N = 3 and D = 2, since this case is particularly easy to visualize. To see why it is easy to visualize, first note that Y is a three-dimensional vector; thus by our discussion above, we can parameterize Y by a vector on the unit circle (i.e. a scalar between 0 and 2π). Second, note that the matrix of left singular vectors U is parameterized by two orthonormal vectors, U·1 and U·2, each on the unit 2-sphere. As both vectors must be orthogonal to 1 ∈ R3, we can parameterize U·1 and U·2 by two orthonormal vectors on the unit circle. We parametrize U·1 by a scalar that determines U·2 up to a rotation of U·2 by π . We fix a rotation for U·2 relative to U·1 so that, for fixed singular values S, the quasiconvexity of L is parameterized by two scalars. Fig. 2 is the resulting visualization. Precisely, to make Fig. 2, we fix an orientation between U·1 and U·2: here, a rotation of π/2 on the unit circle. We create a uniform grid over the unit circles for U·1 and Y . Fig. 2 visualizes the severity of non-quasiconvexity of L as we move over this grid for three different settings of the singular values. To define the severity of non-quasiconvexity, let λworst, λ∗, and λworst-min correspond to the λ maximizing L, the λ minimizing L, and the λ corresponding to the local minimum with largest L, respectively. We then compute
severity := L(λworst-min)− L(λ∗) L(λworst)− L(λ∗)
(3)
If the severity is near 1, then finding the worst local minimum is nearly as bad as selecting the worst possible λ. In Fig. 2 we color pixels according to this measure of severity and do in fact see values near 1. We see that, in general, U , Y , and S have a complicated interaction to determine the quasiconvexity of L. The values of Y that produce quasiconvexity depend on U . No setting of Y or U guarantees quasiconvexity.
Fig. 2 demonstrates the poor performance of a hypothetical optimizer that always finds the worst local minimum. In ??, we show that poor performance is not just hypothetical; in particular, the most commonly used methods for optimizing L – grid search and gradient descent – encounter difficulties due to non-quasiconvexity. Issues due to non-quasiconvexity are also not specific to measuring performance in terms of L nor to optimizing the LOOCV loss. In ??, we show that the K-fold CV loss also suffers from non-quasiconvexity and that computing Eq. (3) with the test loss instead of L also exhibits issues due to non-quasiconvexity.
Despite the complexity of Fig. 2, one trend does seem clear: as the singular values become more similar (moving from left to right in the panels of Fig. 2), the fraction of Y values and U values that correspond to quasiconvexity (dark purple regions) grows. Based on this behavior, one might conjecture that a sufficiently uniform spectrum of the covariate matrix could guarantee the quasiconvexity of L.

heading: 5 Quasiconvexity of L with a nearly uniform spectrum S
text: We now build on the conjecture of the previous section to show that we can, in fact, guarantee quasiconvexity in certain cases. In particular, we will show conditions under which a sufficiently uniform spectrum S of the covariate matrix X guarantees that L is quasiconvex. One might hope that for large N , eventually any Y or U would yield quasiconvexity. However, even when S is exactly uniform, our experiments in Section 6 show that we cannot expect such a statement. Rather, we will devise conditions on Y and U to avoid “extreme” settings for either quantity. With these conditions, our main theorem will show that a sufficiently flat spectrum S does indeed guarantee quasiconvexity of L. When our theorem applies, we can safely terminate an optimization procedure at the first local minimum of L that we discover. We first establish notation and then state our assumptions.
Definition 2. Let θ̂ := ( XTX )−1 XTY be the least-squares estimate. Define the least-squares
residuals Ê := Y −Xθ̂. Let ε̂n be the nth entry of Ê.
Note that θ̂ is well-defined since we have assumed that D < N and that X is full-rank. For the tractability of our theory, all of our assumptions and conclusions will be asymptotic in N ; in particular, our assumptions will use big-O and little-o statements, which are to be taken with respect to N growing large. In our discussion of these assumptions, we assume that D is fixed. There is nothing in our proofs or assumptions that requires D to be fixed; however the validity of our assumptions is not clear if D grows with N , so we do not consider this case here. Since LOOCV is useful precisely for
finite N , we are careful to show in our experiments (Section 6) that these asymptotics take hold for small N .
Our first assumption concerns the magnitude of the residuals Ê. Assumption 1. (1/N) ∑N
n=1 ε̂ 2 n is O(1) (i.e. it does not grow with N ).
This assumption is fairly lax. For example, suppose our linear model is well-specified. In particular, suppose there exists some θ∗ ∈ RD such that yn = ⟨xn, θ∗⟩+ εn where the εn are i.i.d. N (0, σ2) for some σ > 0. Stack the εn into a vector E ∈ RN . Then ∥Ê∥2 = ∥(IN − UUT )E∥2 < ∥E∥2. Since (1/N)∥E∥2 is O(1) with high probability, it follows that Assumption 1 holds with high probability in this well-specified linear model. We emphasize that Assumption 1 depends on the residuals of the least squares estimate, not (directly) on the noise in the observations.
Our next assumption governs the size of the least squares estimate θ̂.
Assumption 2. ∥θ̂∥ is O(1) (i.e. it does not grow with N ).
Again, this is a lax assumption. For example, given any statistical model for which θ̂ is a consistent estimator for some quantity, Assumption 2 holds.
Our next assumption constrains the uniformity of the left-singular value matrix U with rows un ∈ RD. Assumption 3. We have maxn ∥un∥2 := ∥umax∥2 = O(N−p) for some p > 1/2.
Assumption 3 is an assumption about the coherence of the U matrix, a quantity of importance in compressed sensing and matrix completion [Candés and Recht, 2009]. In particular, Assumption 3 requires that the coherence of U decay sufficiently fast as a function of N . Suppose we remove the condition that U have zero-mean columns (see Condition 1 and the discussion after Remark 1) and assume a uniform distribution over valid U (i.e. matrices with orthonormal columns); then Assumption 3 is known to hold with high probability for any p ∈ (1/2, 1) [Candés and Recht, 2009, Lemma 2.2].
There do exist matrices U with orthonormal zero-mean columns that do not satisfy Assumption 3. For instance, take some small N0 (say N0 = 5) and a valid U ′ for this N0. Then, for N > N0, form U by appending 0 ∈ R(N−N0)×D to the bottom of U ′. This construction yields an N ×D matrix U with orthonormal and zero-mean columns for which ∥umax∥2 is constant as N grows. Still, in our experiments in Section 6 and ??, we confirm that, for a uniform distribution over orthonormal U with zero-mean columns, Assumption 3 holds with high probability.
Our final assumption is a technical assumption relating ∥un∥2, ε̂n, and θ̂. Assumption 4. The following quantity is positive and Θ(1) (i.e. is bounded away from zero and does not grow with N ): ∥θ̂∥2 − ∑N n=1 ∥un∥2 ( ⟨un, θ̂⟩2 + 2ε̂2n
) Roughly, this assumptions means that the largest ∥un∥2 and ε̂2n values do not occur for the same values of n. To see this relation, note that Assumption 3 implies ∥θ̂∥2 − ∑ n ∥un∥2⟨un, θ̂⟩2 ≥ (1 − O(N−p))∥θ̂∥2. If we assume that ∥θ̂∥2 = Θ(1) (i.e. Assumption 2 holds and θ̂ does not converge to 0 ∈ RD), then we find ∥θ̂∥2 − ∑ n ∥un∥2⟨un, θ̂⟩2 = Θ(1). So, we need only that∑
n ∥un∥2ε̂2n = o(1) for Assumption 4 to hold; e.g. we need that the largest values of ∥un∥2 and the largest values of ε̂2n typically do not occur for the same values of n.
With our assumptions in hand, we can now state our main theorem. Our theorem relates the uniformity of the spectrum of X to the quasiconvexity of L. As we have shown in Proposition 2, the scaling of the singular values does not matter for the quasiconvexity of L. We therefore take the spectrum to be nearly uniform around 1 ∈ RD. Theorem 1. Consider a series of regression problems with N growing to infinity, where the N th problem uses data (X(N), Y (N)). Assume this sequence satisfies Assumptions 1 to 4. Let the covariate matrix of the N th regression problem have SVD X(N) = U (N) diag(S(N))V (N). There is a N0 > 0 and neighborhood ∆ of 1 ∈ RD such that if N ≥ N0 and the spectrum S(N) ∈ ∆, then L is quasiconvex.
Proof sketch: For one-dimensional functions L, a sufficient condition for quasiconvexity is that for all λ such that L′(λ) = 0, we have L′′(λ) > 0 [Boyd and Vandenberghe, 2009, Chapter 3.4]. We first show L′ can be zero only for a bounded set of λ. We then show that for any λ within this set with L′(λ) = 0, we have L′′(λ) > 0. See ?? for a full proof.
In Section 4, we showed it can be difficult to guess when L is quasiconvex. But Theorem 1 yields one condition that guarantees L is quasiconvex: when X has a nearly uniform spectrum. A natural question then is: when is the spectrum of X nearly uniform? As it happens, a uniform spectrum occurs under standard assumptions, for example, when the xnd are i.i.d. sub-Gaussian random variables.
Definition 3 (e.g. [Vershynin, 2018]). A random variable Q is sub-Gaussian if there exists a constant c > 0 such that E[exp ( Q2/c2 ) ] ≤ 2.
Corollary 1. Take any series of regression problems satisfying Assumptions 3 and 4. Assume the series of regression problems are drawn from a well-specified linear model for some θ∗ ∈ RD: y (N) n = ⟨x(N)n , θ∗⟩ + εn, where εn
i.i.d.∼ N (0, σ2). If σ is sufficiently small, θ̂ is consistent for θ∗, and the entries of the covariate matrices x(N)nd are i.i.d. sub-Gaussian random variables, then L is quasiconvex with probability tending to 1 as N → ∞.
Proof sketch: Assumptions 1 and 2 hold for a well-specified linear model. If the entries of X are i.i.d. sub-Gaussian random variables, standard concentration inequalities imply that its spectrum is nearly uniform with high probability; hence the result of Theorem 1 applies. See ?? for a full proof.

heading: 6 Theorem 1 in practice
text: In Section 5, we established a number of assumptions that we then required in Theorem 1 to prove that L is quasiconvex. A few questions remain about our theorem in practice: (1) how large is the neighborhood ∆, (2) how necessary are our assumptions, and (3) how large do we require N to be? We explicitly answer (1) and (2) with experiments below. (3) is particularly concerning, as regularization has minimal impact when N ≫ D. That is, there is little performance gain by using a regularizer, which removes the need for hyperparameter tuning. To show that our theorem holds when N ∼ D, the majority of our experiments validating Theorem 1 use N that is at most an order of magnitude larger than D.
Throughout our experiments, we check for non-quasiconvexity numerically and use a shortcut formula to compute L that takes advantage of the fact that the right singular vectors of X are V = ID; see ?? for details. The only software dependency for our experiments is NumPy [Harris et al., 2020], which uses the BSD 3-Clause “New” or “Revised” License.
How do we know S is in the neighborhood ∆ in Theorem 1? While our theorem does not give an explicit size of the neighborhood ∆, we can show empirically that ∆ is substantial, even for small to moderate N . We fix D = 5. To generate various spectra of X , we set Sd = eαd/eαD. For α → 0, we get S → 1; we vary α from zero to one to generate spectra of varying distances from uniformity. For each α, we sample 100 left-singular-value matrices U from the uniform distribution over orthonormal U with column means equal to 0; see ?? for how to generate such matrices. We fix a unit-norm θ∗ ∈ RD and for each U , we generate data from a well-specified linear model, yn = ⟨xn, θ∗⟩+ εn, where the εn are drawn i.i.d. from N (0, σ2) with variance σ2 = 0.5. In particular, for each setting of U , we generate 100 vectors Y . For each setting of U and Y , we compute L and check whether it is quasiconvex. In the top left panel of Fig. 3, we report the fraction of problems (out of the 100 ∗ 100 = 10,000 datasets for the corresponding α value) with a non-quasiconvex L versus the distance from uniformity, ∥S − 1∥1. We see that, even for N = 20, the fraction of non-quasiconvex problems quickly hits zero as ∥S − 1∥1 shrinks. To provide a rough practical heuristic, we observe from Fig. 3 that when N = 50 and ∥S−1∥1 ∼ 2, non-quasiconvexity occurs less than 1% of the time. In Fig. 3, ∆ appears to grow slightly with N , so Fig. 3 seems to suggest that we should expect to see little quasiconvexity if ∥S − 1∥1 ≤ 2 and N ≥ 50. In practice, how should we access the spectrum to check this condition? When N and D are small enough we can directly compute S via the singular value decomposition; however, in practice, N or D may be large. If N is large and D is small, we can access the spectrum as the square root of the eigenvalues of XTX . If N is very large, formation of XTX can be expensive; in this case, we suggest the use of randomized sketching to obtain a randomized approximation to the spectrum of X [Woodruff et al., 2014]. Finally, when both N and D are large, we can use spectral density estimation, which gives an estimate of the density of a matrix’s eigenvalues [Lin et al., 2016], and has been shown to successfully scale to large problems [Ghorbani et al., 2019, Yao et al., 2020].
Importance of Assumption 3. We now establish the necessity of Assumption 3 on the decay of ∥umax∥2 with N . To do so, we generate two sets of matrices U as N grows. We generate the first set to satisfy Assumption 3, and we generate the second to violate Assumption 3. In both cases, we will take D = 5 and ten settings of N between N = 10 and N = 300.
To generate the assumption-satisfying matrices U , we proceed as follows. For each N , we draw 500 matrices U from the uniform distribution over orthonormal U matrices with column means equal to 0. For each N , we plot the maximum value of ∥umax∥2 across these 500 U matrices in Fig. 3 (top-right) as a blue dot. We fit a line to these values on a log-log plot, and find the slope is -0.74. This confirms that these matrices satisfy Assumption 3.
To generate assumption-violating matrices U , we proceed as follows. Recall that the smallest N is 10 and D = 5. 100 times, we randomly draw a Usmall ∈ R8×5. We then construct each U by appending N − 8×D zeros to Usmall. For each N , we plot the minimum value of ∥umax∥2 across these 100 U matrices in Fig. 3 (top-right) as an orange dot. Since the minimum of ∥umax∥2 is constant with N , Assumption 3 is violated.
Now we check quasiconvexity. To that end, we randomly select a fixed unit-norm vector θ∗ ∈ RD. For each N , we generate 100 noise vectors E ∈ RN , where the entries En are drawn i.i.d. from N (0, 0.5). For each U and E, we construct Y = USθ∗ +E, where S = 1 ∈ RD. We then compute
the fraction of these (100 ∗ 100 = 10,000) losses L that are non-quasiconvex. In the lower panel of Fig. 3, we plot this fraction against N for both for the assumption-satisfying case (blue) and the assumption-violating case (orange) in Assumption 3. When the assumption is satisfied (blue), we see that the conclusion of Theorem 1 holds: beyond a certain N , there are no settings of U or Y that generate a non-quasiconvex L. We see that in practice, the boundary N is small or moderate (below 50). When the assumption is violated (orange dots), we see that the conclusion of Theorem 1 fails to hold: as N grows, there are still settings of U and Y for which quasiconvexity fails to hold. Finally, we call attention to the vertical axis. Even in the assumption-violating case, the fraction of non-quasiconvex losses is small. It follows that, even for our degenerate U ’s, nearly every combination of noise and U leads to a quasiconvex L. An interesting challenge for future work is to provide a precise characterization of this effect.
Do the ε̂n need to be small (Assumption 1)? Finally, we demonstrate the necessity of Assumption 1, which can be restated as requiring that ∥Ê∥2 grows at most linearly in N . But we also find a suggestion that there may be even more permissive assumptions of interest. To this end, we vary N from 10 to 30. For each N , we generate 4,000 settings of U , each uniform over orthonormal U with column means equal to 0. For each U , we generate 250 unit vectors R such that UTR = 0 (each generated uniformly over such R; see ??). Separately, we consider 60 different norms ν for the vector Ê equally spaced between ν = 0 and ν = 2; these are the same across N . We generate a single unit-norm θ∗ ∈ RD. For each setting of U , R, and ν, we consider the regression problem with covariate matrix U1 and responses Y = U1θ∗ + Ê, where Ê := νR. We record whether L is quasiconvex or not for this problem. For a particular N and particular error-norm ν, we check whether any of the L (across 4,000∗250 = 1,000,000 problems) were non-quasiconvex. Finally, for each N , we find the maximum error-norm νmax,N such that for all ν < νmax,N every regression problem is quasiconvex. We plot a dot at (N, νmax,N ) in the left panel of Fig. 4. We see that in fact the boundary of allowable Ê norms does grow about linearly in N .
Our next plot lends additional insight into how the boundary νmax,N varies with N and is also suggestive of other potential variations on Assumption 1 that might be of interest. In particular, in the right panel of Fig. 4, we consider two particular values of N : N = 10 (blue) and N = 20 (red). For each setting of ν on the horizontal axis, we compute the fraction of non-quasiconvex losses L over all settings of U and R (4,000 ∗ 250 = 1,000,000 problems for each ν). We see that, as expected from the left panel of Fig. 4, the boundary νmax,N is higher for N = 20 than for N = 10. Surprisingly, we also see that at high values of ν, the fraction of non-quasiconvex cases decreases again. We conjecture that in general (i.e. beyond these two particular N ), large amounts of noise leads to little or no non-quasiconvexity. Finally, we note that, as in the bottom panel of Fig. 3, the fraction of non-quasiconvex cases across all ν is low. Again, this small fraction suggests a direction for future work.

heading: 7 Discussion
text: We have shown that the LOOCV loss L for ridge regression can be non-quasiconvex in real-data problems. Local optima need not be global optima. These multiple local optima may pose a practical problem for common hyperparameter tuning methods like gradient-based optimizers, which may get stuck in a local optimum, and grid search, for which upper and lower bounds need to be set.
We proved that the quasiconvexity of L is determined by only a few aspects of a linear regression problem. But we also showed that the quasiconvexity of L is still a complicated function of the remaining quantities, and as of this writing the nature of this function is far from fully understood. Nonetheless, we have provided theory that guarantees at least some useful cases when L is quasiconvex: when the spectrum of the covariate matrix is sufficiently flat, the least-squares fit θ̂ fits the data reasonably well, and the left singular vectors of the covariate matrix are regular. In our experiments, we have confirmed that these assumptions are necessary to some extent: when they are not satisfied, L can be non-quasiconvex. Still, our empirical results make it clear there is more to be explored. We describe some of the directions we believe are most interesting for future work below.
Sharper characterization of when L is quasiconvex. Fig. 2 shows that non-quasiconvexity disappears as the spectrum of X becomes uniform; however, it is clear that there is very regular behavior to the pattern of quasiconvexity even when the singular values of X are non-uniform. We are not able to characterize these patterns at this time but believe these patterns pose a fascinating challenge for future work. Relatedly, our experiments (Section 6) show that when our assumptions are violated, quasiconvexity of L is not guaranteed. However, we have observed that even when L is not guaranteed to be quasiconvex, many settings of U and Y still give quasiconvexity. In many of our experiments, the fraction of non-quasiconvex losses L was extremely small. How many local optima and how bad are they? Without the guarantee of a single, global optimum, it is not clear that we can ever know that we have globally (near-)optimized L. However, notice that our examples in Fig. 1 all have at most two local optima. In simulated experiments, we also typically encountered two local optima in non-quasiconvex losses, although we have not studied this behavior systematically. If L were guaranteed to have only two or some small number of optima, optimization might again be straightforward, even in the case of non-quasiconvexity; an algorithm could search until it finds the requisite number of optima and then report the one with the smallest value of L. Alternatively, one might hope that all local optima have CV loss (and ideally out-of-sample error) close in value to that of the global optimum. Indeed, Kawaguchi [2016] argue that this property holds for certain losses in deep learning. Presumably it is not universally the case that local optima exhibit similar loss since the right panel of Fig. 1 seems to give a counterexample. But it might be widely true, or true under mild conditions. Meanwhile, in the absence of such guarantees, optimization of L should proceed with caution.
Beyond ridge regression. We have shown – in our opinion – surprising non-quasiconvexity for the LOOCV loss for ridge regression. Do similar results hold for simple models outside ridge regression? The regularization parameter in other ℓ2 or ℓ1-regularized generalized linear models is often tuned by minimizing a cross-validation loss. In preliminary experiments, we have found non-quasiconvexity in ℓ2-regularized logistic regression. To what extent do empirical results like those in Fig. 2 or theoretical results like those in Theorem 1 hold for other models and regularizers?

heading: Acknowledgements
text: We thank the anonymous reviewers for suggesting valuable additional experiments and very carefully checking (and correcting) our proofs. WS and TB thank an NSF Career Award and an ONR Early Career Grant for support. MU and ZF gratefully acknowledge support from NSF Award IIS-1943131, the ONR Young Investigator Program, and the Alfred P. Sloan Foundation.

Review summaries:
### Update  I have read the author response; I'd like to thank the authors for the detailed comments and clarifications.  - I am convinced from the authors' response that all of the theoretical issues I raised in my review can be fixed/clarified. Indeed, the response above does exactly this and need only be integrated into the main paper.  - I like the proposed experiments illustrating that GD on the LOO-CV loss can get stuck in local minima. This will be a nice demonstration of why the theoretical results matter in practice --- especially if the resulting test accuracies are quite different between the local min and the global solution.  In light of the response, I have decided to increase my score to a 6.  ==========  The authors conduct a theoretical and empirical investigation of the _quasi-convexity_ of the leave-one-out cross validation (CV) metric for $\ell_2$-regularized least-squares regression (aka ridge-regression). The main theoretical results are: i) counter-examples for general (quasi-)convexity of the CV metric, ii) an enumeration of dataset properties which do not effect quasi-convexity (scaling of the singular values of $X$, scaling of the targets $Y$, and the right singular vectors), and iii) sufficient conditions for quasi-convexity to hold.  The empirical component of the submission validates the sufficient conditions for quasi-convexity --- thus confirming the theoretical developments --- and examines sensitivity of the main result to their violation.  Interestingly, these experiments show that quasi-convexity is highly sensitive to the singular-value spectrum, but relatively robust to irregularity in the residuals and left singular vectors.

This paper demonstrates that leave-one-out cross-validation (LOOCV) risk estimate for ridge regression can have multiple local minima (both in real data and in simulations). It then goes on to prove that under some restrictive set of conditions the LOOCV loss is quasiconvex and only has one local=global minimum.  The paper is easy to follow and fun to read. The fact that LOOCV curve can have local minima was pretty surprising to me (but I may be unaware of some relevant literature). Not sure it's a perfect fit specifically for NeurIPS, but I vote for acceptance.

The manuscript characterizes conditions under which the leave-one-out cross-validation error of ridge regression as a function of the ridge hyper-parameter has quasi-convex shape. The quasi-convexity would make the search for an optimal hyper-parameter value easier due to its single global optimum. Provided that certain mild technical assumptions about the data hold, quasi-convexity tends to emerge when the singular value spectrum of the training data is uniform enough. The effect of the assumptions and their violations are experimentally analysed with simulations.