Bayes beats Cross Validation: Efficient and Accurate Ridge Regression via Expectation Maximization
Authors: Shu Tew
Mario Boley
Daniel F. Schmidt
Venue: NeurIPS 2023 poster (2023)
Keywords: Ridge Regression
Cross validation
Expectation Maximisation
Bayesian methods
TL;DR: We propose a novel EM algorithm for optimizing the ridge regression hyper-parameter that outperforms leave-one-out cross validation in terms of computational complexity and statistical robustness

Abstract: We present a novel method for tuning the regularization hyper-parameter, $\lambda$, of a ridge regression that is faster to compute than leave-one-out cross-validation (LOOCV) while yielding estimates of the regression parameters of equal, or particularly in the setting of sparse covariates, superior quality to those obtained by minimising the LOOCV risk. The LOOCV risk can suffer from multiple and bad local minima for finite $n$ and thus requires the specification of a set of candidate $\lambda$, which can fail to provide good solutions. In contrast, we show that the proposed method is guaranteed to find a unique optimal solution for large enough $n$, under relatively mild conditions, without requiring the specification of any difficult to determine hyper-parameters. This is based on a Bayesian formulation of ridge regression that we prove to have a unimodal posterior for large enough $n$, allowing for  both the optimal $\lambda$ and the regression coefficients to be jointly learned within an iterative expectation maximization (EM) procedure. Importantly, we show that by utilizing an appropriate preprocessing step, a single iteration of the main EM loop can be implemented in $O(\min(n, p))$ operations, for input data with $n$ rows and $p$ columns. In contrast, evaluating a single value of $\lambda$ using fast LOOCV costs $O(n \min(n, p))$ operations when using the same preprocessing. This advantage amounts to an asymptotic improvement of a factor of $l$ for $l$ candidate values for $\lambda$ (in the regime $q, p \in O(\sqrt{n})$ where $q$ is the number of regression targets).

Bayes beats Cross Validation: Fast and Accurate Ridge Regression via Expectation Maximization
Shu Yu Tew
Mario Boley Monash University mario. www.yymonash....ooomoooooooooooe edu
Monash University shu. tew@monash. edu
Daniel E. Schmidt Monash University daniel. dt @monashh edu
Abstract
We present novel method for tuning the regularization hyper-parrane papae--rramtter, A, of ridge regression that is faster to compute than leave-one-out tross-validation (LOOOCV) while yielding estimates of the regression parameters of equal, or par- ticularly in the setting of sparse covariates, superior quality to those obtained by minimising the LOOCV risk. The LOOCV risk can suffer from multiple and bad local minima for finite and thus requires the specification of set of candidate X, which can fail to provide good solutions. In contrast, we show that the proposed method is guaranteed to find unique optimal solution for large enough n, under relatively mild conditions, without requiring the specification of any difficult to determine hyper-parametern parameters. This is based a Bayesian formulation of ridge regression that we prove to have a unimodal posterior for large enough n, allowing for both the optimal À and the regression coefficients to be jointly learned within an iterative expectation maximization (EM) procedure. Importantly, we show that by utilizing loop can be implemented in (minn(1, p)) operations, for input data with n rows and p columns. In contrast, evaluating a single value of À using fast LOOCV costs On min (n, p)) operations when using the same preprocessing. This advantage amounts to an asymptotic improvement factor of / for candidate values for A (in the regime q,p E O(vn) where is the number of regression targets).
appropriate preprocessing step, single iteration of the main EM
Introduction
Ridge regression 25 is one of the most widely used statistical learning algorithms. Given training data X E x=*p and R", ridge regression finds the linear regression coefficients Bx that minimize the fz-reuulrd sum of squared errors, i.c.,
A* = arg min {lly - x,3|| + +||
B
In practice, using ridge regression additionally involves estimating the value for the tuning parameter À that minimizes the expected squared error 33x y) for new data x and y sampled from the same distribution as the training data. This problem is usually approached via the leave-one-out cross-validation (LOOCV) estimator, which can be computed efficiently by exploiting closed-form solution for the leave-one-out test errors for given A. The wide and long-lasting use of the LOOCV approach suggests that it solves the ridge regression problem more or less optimally, both in terms of its statistical performance, as well as its computational complexity.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

EM
1.0
8x
6.0
1e5
2-
0-
1000
500
0.75
0.50
0.25
16a
10a
10" 107
107 30" 102
10.5
10.0
107
107
10" 307
307 20
Figure I: Comparison of LOOCV (with fixed candidate grid of size 100) and EM for setting with sparse covariate vectors of x (x1, ......... (100)) such that r; ~ Ber /00)) i.i.d. and responses y|x~ ~ N(3 x, r²) for increasing noise levels a and sample sizes n. In an initial phase for small n, the number of EM iterations tends to decrease rapidly from initial large number until it reaches a small constant (around 10). In this phase, EM is computationally slightly more expensive than LOOCV (third row) but has a better parameter mean squared error (first row) corresponding less shrinkage (second row). In the subsequent phase, both algorithms have essentially identical parameter estimates but EM outperforms LOOCV in terms of computation by wide margin.
However, in this work, show that LOOCV is outperformed by simple expectation maximization (EM) approach based on a Bayesian formulation of ridge regression. While the two procedures not equivalent, in the sense that they generally do not produce identical parameter values, the EM estimates tend to be of equal quality or, particularly in sparse regimes, superior to the LOOCV estimates (see Figure T). Specifically, the LOOCV risk estimates can suffer from potential multiple and bad local minima when using iterative optimization, or misspecified candidates when performing grid search. In contrast, we show that the EM algorithm finds unique optimal solution for large enough (outside pathological cases) without requiring any hard to specify sapee-paraaeters, which consequence of more general bound on n (Thm. 3.1 that we establish to guarantee the unimodality of the posterior distribution of Bayesian ridge regression-a result with potentially wider applications. In addition, the EM procedure is asymptotically faster than the LOOCV procedure by factor of I where is the number of candidate values for À to be evaluated (in the regime p. E ((yn) where p. 9. and are the number of covariates, target variables, and data points, respectively). In practice, even in the usual case of g LOOCV computationally by an order of magnitude as we demonstrate on a test suite datasets from the UCI machine learning repository and the UCR time series classification archive.
and / O(1), the EM algorithm tends to outperform
While the EM procedure discussed in this paper is based on a recently published procedure for learning sparse linear regression models (44), the adaption of this procedure to ridge regression has not been previously discussed the literature. Furthermore, direct adoption would lead to main loop complexity of O(²²) that is uncompetitive with LOOCV. Therefore, in addition to evaluating the empirical accuracy and efficiency of the EM algorithm for ridge regression, the main technical contributions of this work are to show how certain key quantities can be efficiently computed from either a singular value decomposition decomposition of the Gram matrix X "X, when n > This results in an E-step of the algorithm in time (r) where r = min(n, p), and an M-step found closed form and solved time O(1), yielding altra-fast main loop for the EM algorithm.
the design matrix, when p 2 n, or an eigenvalue

---

Table Time complexities of algorithms: m max(n, p). r min(n, for LOOCV, k number of EM iterations and is the number of the target variables.
(n, p). number of candidate
METHOD
MAIN LOOP
PRE-PROCESSING
OVERALL (p.9 E O(vn)
NAIVE ADAPTION OF EM
O(p''g)
O(p2n)
O(kn²) Okn + n²) O(Ir²)
PROPOSED BAYESEM
(krrg) (lverg)
(
FAST LOOCV
O(ar*))
These computational advantages result in an algorithm that is computationally superior efficient, optimized implementations of the fast LOOCV algorithm. Our particular implementation of LOOCV actually outperforms the implementation in scikit-leard by approximately a factor of two by utilizing similar preprocessing to the EM approach. This enables an O(nr) evaluation for single (which is still slower than the (r)) evaluation for our new EM algorithm; see Table for an overview of asymptotic complexities of LOOCV and EM), and may be of interest to readers by itself. Our implementation of both algorithms, along with all experiment code, are publicly available in the standard package ecosystems of the R and Python platforms, as well as on GitHut
In the remainder of this paper, we first briefly survey the literature of ridge regression with an emphasis on the use of cross validation (Sec. - Based on the Bayesian interpretation of ridge regression, we then introduce the EM algorithm and discuss its convergence (Sec. B.. Finally, develop fast implementations of both the EM algorithm and LOOCV (Sec. and compare them empirically (Sec.5) 55)
2
Ridge Regression and Cross Validation
(also known as f--regularization) is a popular method for estimation and prediction in linear models. The ridge regression estimates are the solutions the penalized least- squares problem given in (I). The solution to this optimization problem is given by:
Ridge regression 25
B = (x + AI,))
(2)
When À -> O, the ridge estimates coincide with the minimum & norm least squares solution 22 31 which simplifies to the usual least squares estimator in cases where the design matrix X has full column rank (i.e. X is invertible). Conversely, as X + oc, the amount of shrinkage induced by the penalty increases, with the resulting ridge estimates becoming smaller for larger values of X. Under fairly general assumptions [26), including misspecified models and random covariates of growing dimension, the ridge estimator is consistent and enjoys finite sample risk bounds for all fixed 2 0, i..., it converges almost surely to the prediction risk minimizer, and its squared deviation from this minimizer is bounded for finite n with high probability. However, its performance can still vary greatly with the choice of X; hence, there is need estimate the optimal value from the given training data.
Earlier approaches to this problem [e.g. 74 9 23
34 48 rely explicit estimate 30
30
of the (assumed homoskedastic) noise variance, following the original idea of Hoerl and Kennard 25 However, estimating the noise variance can be problematic, especially when p is not much smaller than n [8] 20, 24 (46). More recent approaches adopt model selection criteria to select the
optimal À without requiring prior knowledge or estimation of the noise variance. These methods involve minimizing selection criterion of choice, such as the Akaike information criterion (AIC) (L] Bayesian information criterion (BIC) (40), Mallow's conceptual prediction (C,,) criterion most commonly, cross validation (CV)
[3], and,
A particularly attractive variant of CV is leave-one-out cross validation (LOOCV), also referred to as the prediction error sum of squares (PRESS) statistic the statistics literature
CV (A)
(A)
(JA 04
n
3

---

where ú = = and 3 denotes the solution to ( when the i-th data point (x,, 34) is omitted. LOOCV offers several advantages over alternatives such as 10-fold CV: deterministic, nearly unbiased 47.. and there exists an efficient "shortcut" formula for the LOOCV ridge estimate [36
R
(A) = Σ (a())
C
where (A) = X(XXX + AI is the regularized "hat", or projection, matrix and = y H(Jyy the residuals of the ridge fit using all n data points. As it only requires the diagonal entries the hat matrix, Eq. (4) allows for the computation of the PRESS statistic with the same time complexity O(p² + np"
single ridge regression fit.
Moreover, unless p/n + 1, the LOOCV ridge regression risk as function of À converges uniformly (almost surely) to the true risk function on [0, oo) and therefore optimizing it consistently estimates the optimal À 22 36 However, for finite n, the LOOCV risk can be multimodal and, even worse, there exist local minima that are almost as bad as the worst À [42]. Therefore, iterative algorithms like gradient descent cannot be reliably used for the optimization, giving theoretical justification for the pre-dominant approach of optimizing over finite grid of candidates L (A1,... A²). Unfortunately. despite the true risk function being smooth and unimodal, a naïvely chosen finite grid cannot be guaranteed contain any candidate with risk value close to the optimum. While this might not pose a problem for small n when the error estimating the true risk via LOOCV is likely large, it can potentially become a dominating source of error for growing and p. Therefore, letting ! grow moderately with the sample size appears necessary, turning it into a relevant factor in the asymptotic time complexity.
As a further disadvantage, LOOCV (or CV in general) is sensitive to sparse covariates, as illustrated in Figure where the performance of LOOCV, relative to the proposed EM algorithm,, degrades as the noise variance o² grows. In the sparse covariate setting, a situation common in genomics, the information about each coefficient concentrated in only few observations. As LOOCV drops observation to estimate future prediction error, the variance of the CV score can be very large when the predictor matrix is very sparse, the estimates depend on only small number of the remaining observations. In the most extreme case, known as the multiple means problem [27]. X In. and all the information about each coefficient concentrated in single observation. In this setting, the LOOCV score reduces to Σ y². and provides no information about how to select A. In contrast, the proposed EM approach explicitly ties together the coefficients via the probabilistic Bayesian interpretation of À as the inverse-variange of the unknown coefficient vector. This "borrowing of strength" means that the procedure provides sensible estimate of À even in the case of multiple means (see Appendix A).
3 Bayesian Ridge Regression
The ridge estimator (2) has well-known Bayesian interpretation; specifically, if we assume that the coefficients are a priori normally distributed with mean zero and common variance 7²0² we obtain a Bayesian version of the usual ridge regression procedure, i.c..
y 0, Nn (XB3, d²In),
3 3172.0 02 Np (0, I²O2I,)) 0² do
(5)
dd?
where x--) is an appropriate prior distribution assigned to the variance hyperparameter 72. For
given 7 0 and a the conditional posterior distribution of (3 is also normal 332
B B|
0 y Np A²
B, A
(6)
A, (X +7
+ 21,,
where the posterior mode (and mean) 3, is equivalent to the ridge estimate with penalty X = 1/72 (we rely the variable name in the notation 3, to indicate whether it refers to ( or 2)).

---

To estimate the 7 hyperparameter in the Bayesian framework, we first must choose a prior distribution for the hypervariance r². We assume that no strong prior knowledge on the degree of shrinkage of the regression coefficients is available, and instead assign the recommended default beta-prime prior distribution for 7² [15 38 with probability density function:
Shrinkage Prior
π (T² (T²) 0 (1+ +1++))-a=b
a a> 2,0b 0,
B(a,b)
where B(a,b) b) is the beta function. Specifically, we choose = b 1/2, which corresponds to standard half-Cauchy prior on т. The half-Cauchy a heavy-tailed, weakly informative prior that is frequently recommended as a default choice for scale-type hyperparameters such as T [38]. Further, this estimator very insensitive to the choice of a or b. As demonstrated by Theorem 6.1 in (6. the marginal prior density over B. (())))))(()) d) = x (33a,, b) has polynomial tails in [|3]| for all a > 0, > 0, and has Cauchy or heavier tails for b S 1/2. This type of polynomial-taile prior distribution over the norm of the coefficients is insensitive to the overall scale of the coefficients, which is likely unknown a priori. This robustness is in contrast to other standard choices of prior distributions for 7² such as the wverse-gamme distribution [e.g.. 35 which are highly sensitive to the choice of hyperparameters (38.
The asymptotic properties of the posterior distributions in Gaussian linear models (5) have been extensively researched [8] [6] I7 45). These studies reveal that in linear models, the posterior distribution of /3 is consistent, and converges asymptotically to normal distribution centered on the true parameter value. When is fixed, this assertion can be established through the Vernstein-Von Mises theorem 45 Sec. 10.2]. Our specifie problem (5) satisfies the conditions for this theorem to hold: 1) both the Gaussian-linear model piy|3, o²) and the marginal distribution f p(y],9, a² )= (3)+2)d3
Unimodality and Consistency
p are identifiable; 2) they have well defined Fisher information matrices; and 3) the priors over 33 and T are absolutely continuous. Further, these asymptotic properties remain valid when the number of predictors pn, allowed to grow with the sample size n at a sufficiently slower rate [6].
The following theorem (see the proof in Appendix B) provides simple bound the number of samples required to guarantee that the posterior distribution for the Bayesian ridge regression hierarchy given by (5) has only one mode outside small environment around zero
Theorem 3.1. Let > O, and let "/n be the smallest eigenvalue of xT /n. 1f > 0 and e > 4/(7))) then the joint posterior p(3,o2, ,e,,,11,11 ,?ly) has a unique mode with 7² 2 e. In particular if In 2 en aa with a < 1 and e > then there is a unique mode with 72 2 e if > (4/(ce)) 1(1 -a).
In other words, all sub-optimal non-zero posterior modes vanish for large enough if the smallest eigenvalue of X grows at least proportionally to some positive power of n. This is a very mild assumption that is typically satisfied in fixed as well random design settings, e.g., with high probability when the smallest marginal covariate variance is bounded away from zero
Expectation Maximization Given the restricted unimodality of the joint posterior 5 for large enough n. in conjunction with its asymptotic concentration around the optimal Bo. estimating the model parameters via an EM algorithm appears attractive, as they are guaranteed converge to an exact posterior mode. In particular, in the non-degenerate case that Go a 0, there exist 72 = c² > O. such that for large enough, but finite n, the posterior concentrates around (3,,r?), and thus Bg is identified by EM if initialized with a large enough ?
Specifically,, we use the novel approach [44 in which the coefficients B are treated as "missing data", and 72 and o² as parameters to be estimated. Given the hierarchy ( the resulting Bayesian EM algorithm then solves for the posterior mode estimates of 3 by repeatedly iterating through the following two steps until convergence:
E-step. Find the parameters of the Q-finction, i.c., the expected complete negative log- posterior (with respect to B3, conditional on the current estimates of f? and o?. and the observed data y:
2011..0?)) of)) = Ep[- -loogp(B,,
o2 Iy) /7,
x7.x]
n1 /n+p+2
log 0² ESS
p + 1
ESN
20 +
log 2 +
20 + log(1 + +73) (8)
2
5

---

where the quantities to be computed are the (conditionally) expected sum of squared errors ESS
E [ll x f.. of] and the expected squared norm ESN E [||3]] 2 . o.]. Denoting by
tr(-) the trace operator, one can show (see Appendix C) that these quantities can be computed as
ESS = Ily X 3-11
and ESN = (rtr(A, + 1|
M-step. Update the parameter estimates by minimizing the Q-function with respect to the shrinkage hyperparameter 7² and noise variance x² i.e.,
f = arg arg min {Q (F2 1 02:
(10)
,9,0²
Instead of numerically optimizing the two-dimensional Q-function (10). we can derive closed-form solutions for both parameters by first finding 6 (T()), i.e., the update for o² then substituting this into the Q-function. This yields a Q-function that is no longer dependent on o², and solving for 72 is straightforward. The resulting parameter updates in the M-step given by:
function of 72, and
0² ESSS + ESN
(n - 1)ESN - (1 + p) p)ESS + Vg and ²²=
(11)
(n p + 2)
(6 + 2p)
where g (4n + ESN (3 + ESS + ((1 ESN + (p - 1)ESS)² The derivations of these formulae are presented in Appendix D.
From [II. we see that updating the parameter estimates in the M-step requires only constant time. Therefore, the overall efficiency of the EM algorithm is determined by the computational complexity of the E-step. Computing the parameters of the Q-functions directly via (9) requires inverting A,, resulting in O(p²) operations. In the next section, we show how to substantially improve this approach via singular value decomposition.
Fast Implementations via Singular Value Decomposition
To obtain efficient implementations of the E-Step of the EM algorithm as well as of the LOOCV shorteut formula, one can exploit the fact that the ridge solution preserved under orthogonal transformations. Specifically, let = min(n, p) and m max p) and let UEVT = X be compact singular value decomposition (SVD) of X. That is, U E Raxr and V E Rxxr are semi- orthonormal column matrices, i.e., U In and VTV = Ip. and Σ = diag(s), s,) E grxr a diagonal matrix that contains the non-zero singular values s = (s1..... 8,) of X. With this decomposition, and an additional (nr) pre-processing step to compute the ridge solution a E IR" for given T with respect to the rotated inputs XV in time O(r) via
EU'y, y, we can compute
å = -'"""U"' + EuTy- (E² + - = (1/(s) + - )) Oc (12)
where a @ b denotes the element-wise Hadamard product of vectors and b. The compact SVD itself can be obtained in time O(mr²) via an eigendecomposition of either X = VV²VT in case 2 p XXT = UT in case n < p followed by the computation of the missing U = XVE- or = -U
In summary, after an O(mr²) pre-processing step, we can obtain rotated ridge solutions for an individual candidate 7 in time O(r). Moreover, for the optimal 7*, we can find the ridge solution 3, = Vå with respect to the original input matrix via an Oprr) post-processing step. Below show how the key statisties that have to be computed per candidate T (and σ)) can be computed efficiently based on o,, the pre-computed c, and SVD. For the EM algorithm, these are the posterior squared norm and sum of squared errors, and for the LOOCV algorithm, this is the PRESS statistic. While the main focus of this work is the EM algorithm, the fast computation of the PRESS shortcut formula appears to be not widely known (e.g., the current implementation in both scikit-learn and gl mnet t do not use it) and may therefore be of independent interest.
ESN For the posterior expected squared norm ESNN o² (rrr(A,))) + ||,,11 we first observe that || 1| = | = || and then note that the trace can be computed as
tr( (A,AA)) = tr V p E + T -VV!) p
72 max(p n, 0) Σ 1/(5/+ +T
(13)
j=1

---

where in the first equation we denote by Vp. Dn the full matrices of eigenvectors and eigenvalues of X (including potential zeros), and in the second equation used the eyclical property of the trace. Thus, all quantities required for ESN can be computed in time ((r) given the SVD and &
ESS For the posterior expected sum of squared errors ESS = - 1 x
we can compute the residual sum of squares term via
lly x = |yy 2y UU + |UE 11
= | - - || à.1??
(14)
where use B, Vá and X = UEVT in the first equation and the orthonormality of U and the definition of c EUT y in the second. Finally, for the trace term, we find that
ur(x T X = (r(VE'VT(V(E" + II, VT)-)
Σ +
(15)
PRESS The shortcut formula of the PRESS statistic 4 for candidate À requires the computation
of the diagonal elements of the hat matrix H{) = X
X + AIp) "XX" as well as the residual
vector y Hy. With the SVD, the first simplifies to
H(A) UΣ(Σ'
+ XL,) EU"
s
U diag
UT
+
where we use the fact that diagonal matrices commute. This allows to compute the desired diagonal elements hi in time (r) via
(A) = Σ /553 +
(16)
1=1
where May denotes the elements of U. Computing the residual vector is easily done via the rotated ridge solution e = y UΣ However, this still requires O(nr) operations, simply because there are n residuals to compute.
Thus, in summary, by combining the pre-processing with the fast computation of the PRESS statistic, we obtain an overall O(mr² +(qur) implementation of ridge regression via LOOCV where denotes the number of candidate À and q the number of regression target variables. In contrast, for the EM algorithm, by combining the fast computation of ESS and ESN, we end up with an overall complexity of Omrr kgr) where & denotes the number of EM iterations. If we further assume that A = o(n). which is supported by experimental results, see Sec. 5) and that both g.p = O(yn) there is an asymptotic advantage of factor of 1 of the EM approach. This regime is common in settings where more data allows for more fine-gran fine-grained input as well as output measurements, e.g., in satellite time series classification via multiple target regression [I 37). All time complexities are summarized in Tab. [ and detailed pseudocode for both the fast EM algorithm and the fast LOOCV algorithm is provided in the Appendix (see Table 3 and 4).
Empirical Evaluation
In this section, we compare the predictive performance and computational cost of LOOCV against the proposed EM method. We present numerical results on both synthetic and real-world datasets. To implement the LOOCV estimator, we use predefined grid, L (Ai...., .......,, A1). We use the two most common methods for this task: (i) fixed grid arbitrarily selecting a very small value as Amin- large value Amaxx and construct sequence of values from Amax to Amin on log scale; (ii) data-driven grid - find the smallest value of Amax that sets all the regression coefficient vector to zero 2 (i.e. _ 0), multiply this value by a ratio such that Amin kdmax and create sequence from
Forr ridge regression, Amax = cc. Following the glmnet package, the sequence of À is derived for 0.001. The penalty function used by gimnet A((1 - -a]]|[13 [:1]B1]1). where = 0 corresponds to ridge regression.
7

---

p=
p=100
p==400
p=
10"1
EM
Cv fx
CV.slm
3 10-4
107
6000
4.00
2000
1.4
0.5
0.0
500
1000
15.00
1000
1500
1500
Figure 2: Comparison of EM to LOOCV variants for increasing n and p for settings with N, Σ) and y|x N(x B, 0.25) with random 3 ~ N, I) and E W,(T,)).
Amax Amia on log scale. The latter method is implemented in the glmnet package in combination with an adaptive K coefficient
0.0001
if 2
0.01
otherwise
which we replicate here as input fast LOOCV algorithm (Appendix, Table 4) to efficiently
recover the glmnet LOOCV ridge estimate. 3
We consider fixed grid of A (10-10 1010) and the grid based the glmnet heuristic; in both cases, we use a sequence of length 100. The latter is a data-driven grid, so we will have different penalty grid for each simulated real data set. Our EM algorithm does not require predefined penalty grid, but it needs a convergence threshold which we set to be E = 10 All experiments in this section are performed in Python and the R statistical platform. Datasets and code for the experimental results is publicly available. As is standard in penalized regression, and without standardized the data before model fitting. This means that the predictors
any loss of generality, are standardized to have zero mean, standard deviation of one, and the target has a mean of zero, i.e.,
the intercept estimate is simply Bu-- = (1/n) Σ y1.
5.1 Simulated Data
In this section, use simulation study to investigate the behavior of EM and LOOCV function of the sample size, n, and two other parameters of interest: the number of covariates and the noise level of the target variable. In particular, we are interested in the parameter estimation performance, the corresponding X-values, and the computational cost. To gain further insights into the latter, the number of iterations performed by the EM algorithm is of particular interest, we do not have quantitative bounds for its behavior. We consider two settings that vary in the level of sparsity and correlation structure of the covariates. The first setting (Fig DD assumes i.i.d Bernoulli distributed covariates with small success probabilities that result in sparse covariate matrices, while the second setting (Fig. 2 assumes normally distributed covariates with random non-zero covariances. both cases, the target variable is conditionally normal with mean x" B for random 3 drawn from standard multivariate normal distribution.
I Looking at the results, a common feature of both settings is that the computational complexity of the EM algorithm is non-monotone function in n. In contrast LOOCV, the behavior of EM shows distinctive phases where the complexity temporarily decreases with n before it settles into the, usually expected, monotonically increasing phase. As can be seen, this due to the behavior of the number of iterations k which peaks for small values of n before it declines rapidly to small
"gmmne LOOCV is computed directly by model refitting via coordinate- wisee descent which can be slow.

---

constant (around 10) when the cost of the pre-processing begins to dominate. The occurrence of these phases is more pronounced for both growing p and growing o. This behavior is likely due to the convergence to normality of the posterior distribution as the sample size n oo, with convergence being slower for large p-
An interesting observation is that CV with the employed gimnet grid heuristic fails, in the sense that the resulting ridge estimator does not appear to be consistent for large p in Setting This due to the minimum value of À produced by the glmnet heuristic being too large, and the resulting ridge estimates being overshrunk. This clearly underlines the difficulty of choosing a robust dynamic grid- a problem that our EM algorithm avoids completely.
5.2 Real Data
We evaluated our EM method on 24 real-world datasets. This includes 21 datasets from the UCI machine learning repository [5] (unless referenced otherwise) for normal linear regression tasks and 3 time-series datasets from the UCR repository [10] for multitarget regression tasks. The latter is a multilabel classification problem in which the feature matrix was generated by the of-te-of-thh-ar HYDRA [2] time series classification procedure (which by default uses LOOCV ridge regression for classification), and we train q ridge regression models in a one ver-usrus--ll fashion, where q is the number of target classes. The datasets were chosen such that they covered a wide range of sample sizes, n, and number of predictors, p. We compared our EM algorithm against the fast LOOCV in terms of predictive performance, measured in R2 (and classification accuracy) on the test data, and computational efficiency.
Our linear regression experiments involve settings: (i) standard linear regression; (ii) second- order multivariate polynomial regression with added interactions and second-degree polynomial transformations of variables, and (iii) third-order multivariate polynomial regression with added three-way interactions and cubic polynomial transformations. For each experiment, repeated the process 100 times and used a random 70/30 train-test split. Due to memory limitations, we limit our design matrix size to a maximum of 35 million entries. If the number of transformed predictors exceeded this limit, we uniformly sub-sampled the interaction variables to ensure that pp I 35000000/ (0.7n., and then fit the model using the sampled variables. Note that always keep the original variables (main effects) and sub-sampled the interactions. In the case of multitarget regression, performed random 70/30 train-test split and repeated the experiment 30 times. To ensure efficient reproducibility of our experiments, we set a maximum runtime of 3 hours for each dataset. Any settings that exceeded this time limit were consequently excluded from the result table.
Table 2 details the results of our experiments; specifically, the ratio of time taken to run fast LOOCV divided by the time taken to run our EM procedure (2), and the R² values obtained by both methods on the withheld test set. The number of features, p. and observations, n recorded are values after data preprocessing (missing observations removed, one-hot encoding transformation, etc.). The results demonstrate that our EM algorithm can be up to 49 times faster than the fast LOOCV, with the speed-ups becoming more apparent the sample size n and the number of target variables q increases. In addition, we see that this advantage in speed does not come at a cost in predictive performance, as our EM approach is comparable to, if not better than, LOOCC in almost all cases (also see Appendix, Figure 1, in which most of I²2 values are distributed along the diagonal line).
An interesting observation is that LOOCV using the fixed grid can occasionally perform extremely poorly (as indicated by large negative R² values) while LOOCV using the glmnet grid does not seem to exhibit this behavior. This appears likely to be due to the grid chosen using the glmnet heuristic. Its performance is artificially improved because it is unable to evaluate sufficiently small values of À and is not actually selecting the very small À value that minimizes the LOOCV score. The incorrectly large À values are providing protection in these examples from undershrinkage.
Conclusion
The introduced EM algorithm is a robust and computationally fast alternative to LOOCV for ridge regression. The unimodality of the posterior guarantees a robust behavior for finite n under mild conditions relative to LOOCV grid search, and the SVD preprocessing enables an overall faster computation with an ultra-fast O(k p)) main loop. Combining this with procedure such as orthogonal lqaat-s to provide a highly efficient forward selection procedure is a promising avenue for future research. As the fuuuciio is expected negative logttosterior, it offers a score

---

Table 2: Real datasets experiment results. The first column is the dataset (abbreviated, refer to Appendix E for the full name); the number of target variables q, the training sample size the raw number of features p; T the ratio of time tcv hcv/tea ttM p is the number of features including interactions; EM, Fix, and GLM are the R values the test data for the three procedures.
LINEAR
2ND ORDER
3RD ORDER
DATASET
p
T
EM
GLM
T
EM
GLM
p
T
EM
GLM
TWITTER
408275
BLOG
39355
275
13
0.46
0.46
0.46
804
9.1
0.51
0.51
0.51
CT SLICES
37450
379
12
0.86
0.86
0.86
930
7.7
0.92
0.91
0.92
TOMSHW
19725
96
0.63
0.63
0.63
1775
6.5
0.71
0.71
0.71
NPD COM
8353
13
13
0.84
0.84
0.84
104
15
1.00
1.00
1.00
559
8.6
1.00
1.00
1.00
NPD TUR
PT MOTOR
4112
19
13
0.15
0.15
0.15
208
12
0.25
0.19
0.21
1539
3.9
-1.09
0.01
0.04
TOTAL
19
1539
ABALONE
2923
9
13
0.53
0.53
0.53
51
16
0.38
0.35
0.50
209
12
0.28
0.12
0.12
CRIME
1395
99
14
0.66
0.66
0,66
5049
1.3
0.67
-0.74
0,66
17652
1.1
0.66
-0.22
0.60
AIRFOIL
1052
0.51
0.51
0.51
0.62
0.62
0.73
0.73
0.73
STUDENT
730
39
12
0.18
0.18
0.18
801
3.8
0.19
-0.89
0.16
10693
1.1
0.19
-6.22
-0.18
F.
361
12
2.4
-0.01
-0.01
-0.03
295
0.5
-0.01
-0.01
-0.14
1984
0.3
-0.01
-50.6
-0.45
H.USSI
13
0.71
0.71
0.84
2.2
0.83
-362
0.83
FACEBOOK
346
15
15
0.91
0.91
0.91
167
6.6
-5.09
-26.4
-3,99
1087
1.9
-2.53
-2B4
-5.76
DIABETES
309
10
13
0.49
0.49
0.49
65
6.5
0.49
0.48
0.48
285
2.1
0.47
0.47
0.47
R.ESTATE
6
13
0.56
0.56
0.56
0.65
0.65
5.5
0.65
0.65
0.65
A.MPG
278
8
0.81
0.81
0.81
8.3
0,86
0.86
0.86
119
4.1
0,86
0.86
0.86
YACHT
215
7
16
0.97
0.97
0.97
11
0.97
0.97
0.97
6.1
0.98
0.98
0.98
A.MOBILE
111
25
6.1
0.90
0.89
0.89
1076
1.7
0.90
0.89
12924
0.5
0.88
-184
0.83
EYE
84
200
2.7
0.50
0.26
0.45
20300
1.3
0.19
0.19
0.19
RIBO
49
4088
2.3
0.64
0.64
0.64
CROP (24)
11760
3072
FLECD (7)
11645
4096
9.2
0.88
0.88
0.89
STARL (3)
6465
7168
2.1
0.98
0.60
0.98
This dataset is not from the UCI repository, Data references can be found the appendix. TTme--sriee dataset from the UCR repository. EM, Fix, and GLM the classilication accuracy the data
on which the usefulness of predictors themselves may be assessed, i.e., model selection criteria, resulting in a potentially very accurate and fast procedure for sparse model learning. An important open problem the theoretical analysis of the expected number of EM iterations that required for convergence. The empirical evidence suggests that & converges to a constant, and is thus negligible in the asymptotic time complexity. This in alignment with the convergence of the posterior multivariate normal distribution. However, such intuitive and empirical arguments cannot replace a rigorous worst-case analysis.
Acknowledgments and Disclosure of Funding
We thank the anonymous reviewers for their valuable feedback that led a substantially improved theoretical analysis. This work was supported by the Australian Research Council (DP210100045),
References
[1] Hirotugu Akaike. new look at the statistical model identification. Selected Papers of Hirotugu Akaike, pages 215-222. Springer, 1974.
[2] Mahdi Alkhamisi, Ghadban Khalaf, and Ghazi Shukur. Some modifications for choosing ridge parameters. Communications in Statistics and Methods, 35(11):2005-2020, 2006.
[3] David M Allen. The relationship between variable selection and data agumentation and a method for prediction. technometrics, 16(1):125-112 1974.
[4] Sylvain Arlot and Alain Celisse. A survey of cross-validation procedures for model selection. 2010.
[5] Arthur Asuncion and David Newman. UCI machine learning repository, 2007.
10

---

[6] Ole Crindorff-Nielsene John Kent, and Michael Sørensen. Normal variance mean mixtures and distributions. International Statistical Review/Revue Internationale de Statistique, pages 145-159, 1982.
[7] Satish Bhat and Vidya Raju. class of generalized ridge estimators. Communications in Statistics Simulation and Computation, 46(7):5105-5112, 2017. doi: 10. 10.100000310011 0.260. 1144765, URL https: : http:://..i.........eet//o.o 201//10.10800066009122 2016 1144765
[8] Dominique Bontemps. Bernstein-von Mises theorems for Gaussian regression with increasing number regressors. 2011.
[9] Peter Bühlmann, Markus Kalisch, and Lukas Meier. High-dimensional statisties with view toward applications in biology. Annual Review of Statistics and Its Application, 1:255-278, 2014.
[10] Hoang Anh Dau, Anthony Bagnall, Kaveh Kamgar, Chin-Chia Michael Yeh, Yan Zhu, Shaghayegh Gharghabi, Chotirat Ann Ratanamahatana, and Eamonn Keogh. The UCR time series archive. IEEE/CAA Journal of Automatica Sinica, 6(6):1293-1305, 2019.
[11] Angus Dempster, François Petitjean, and Geoffrey I Webb. ROCKET: Exceptionally fast and accurate time series classification using random convolutional kernels. Data Mining and Knowledge Discovery, 3455)1445--145, 2020.
[12] Angus Dempster, Daniel F Schmidt, and Geoffrey Webb. HYDRA: kernels for fast and accurate time series classification. Data Mining and Knowledge Discovery, 2023.
Competing convolutional
[13] Bradley Efron, Trevor Hastie, lain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of statistics, 32(2):407-499 2004.
[14] AK Md Ehsanes Saleh and BM Golam Kibria. Performance of some new preliminary test ridge regression estimators and their properties. Communicutions in Statistics-Theory and Methods, 22(10):2747-2764 1993.
[15] Andrew Gelman. Prior distributions for variance parameters in hierarchical models. Buyesian Analysis, 1(3):515-533, 2006.
[16] Subhashis Ghosal. Asymptotic normality of posterior distributions in high- high-dimensional linear models. Bernoulli, pages 315-331, 1999.
[17] Subhashis Ghosal, Jayanta K Ghosh, and Aad W Van Der Vaart. Convergence rates of posterior 2000.
distributions. Annals of Statistics, pages 500-531, 2000.
[18] Gene H Golub, Michael Heath, and Grace Wahba. Generalized cross cross- validation as a method for choosing good ridge parameter. Technometrics, 21(2):115-223 1979.
[19] Ramadan Hamed, Ali EL Hefnawy, and Aya Farag. Selection of the ridge parameter using mathematical programming. Communications in Statistics-Simulatinn and Computation, 42(6): 1409-1432, 2013.
[20] Richard J Hanson. A numerical method for solving Fredholm integral equations of the first kind using singular values. SIAM Journal on Numerical Analysis, 8(3):616 6222,, 1971.
[21] David Harrison Jr and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1):81-00,, 1978.
[22] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani. Surprises in high- dimensional ridgeless least squares interpolation. Annals of statistics, 50(2):949, 2022.
[23] Ali El Hefnawy and Aya Farag. A combined nonlinear programming model and Kibria method for choosing ridge parameter regression. Communications Statistics-Simul Simulation and Computation, 43(6): 1442-1470, 2014.

---

[24] John W Hilgers. On the equivalence of regularization and certain reproducing kernel Hilbert space approaches for solving first kind problems. SIAM Journal on Numerical Analysis, 13(2): 172-184, 1976.
[25] Arthur E Hoerl and Robert W Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55 67, 1970.
[26] Daniel Hsu, Sham M. Kakade, and Tong Zhang. Random design analysis of ridge regression. In Shie Mannor, Nathan Srebro, and Robert C. Williamson, editors, Proceedings of the 25th Annual Conference on Learning Theory, volume 23 of Proceedings of Machine Learning Research, pages 9. 9.,4,,4, Edinburgh, Scotland, 25-27 Jun 2012. PMLR.
[27] W JAMES. Estimation with quadratic loss. In Proc. 4th Berkeley Symp. on Math. Statist, and Prob., 1961, 1961.
[28] Lawless JF, A simulation study of ridge and other regression estimators. Communications in Statistics-therry and Methods, 5(4):307-323, 1976.
[29] Ghadban Khalaf and Ghazi Shukur. Choosing ridge parameter for regression problems. 2005.
[30] BM Golam Kibria. Performance of some new ridge regression estimators. Communications in PaatisticsSSmullatie and Computation, 32(2):419 -35.455,999 2003.
[31] Dmitry Kobak, Jonathan Lomond, and Benoit Sanchez. The optimal ridge penalty for worldworld high-dimensional data can be zero or negative due to the implicit ridge regularization. The Journal of Machine Learning Research, (1()):6663-68 6678,, 2020.
[32] D. V. Lindley and A. M. Smith. Bayes estimates for the linear model. Journal of the Royal Statistical Society (Series B), 34(1):1-41, 1972.
[33] Colin L Mallows. Some comments on Cp. Technometrics, 42(1):87-94 2000.
[34] Gisela Muniz and BM Golam Kibria. On some ridge regression estimators: An empirical comparisons. Communications in
and Computation@.
(tttaaatiooon333)666111111 38(3):621-63))
2009.
[35] Trevor Park and George Casella. The Bayesian lasso. Journal of the American Statistical Association, 103(482):681 1-66, 2008
[36] Pratik Patil, Yuting Wei, Alessandro Rinaldo, and Ryan Tibshirani. Uniform consistency of cross-validation estimators for nigh-dimensional ridge regression. In International Conference on Artificial Intelligence and Staristies, pages 3178-3186, PMLR, 2021.
[37] François Petitjean, Jordi Inglada, and Pierre Gançarski. Satellite image time series analysis under time warping. IEEE transactions on geoscience and remote sensing, 50(8):3081 11.301110995,,, 2012.
[38] Nicholas G Polson and James G Scott. On the half-Cauchy prior for global scale parameter. Bayesian Analysis, 7(4)8887- 2012.
[39] Todd E Scheetz, Kwang- Youn A Kim, Ruth E Swiderski, Alisdair R Philp, Terry A Braun, Kevin L Knudtson, Anne M Dorrance, Gerald F DiBona, Jian Huang, Thomas L Casavant, et al. Regulation of gene expression in the mammalian eye and its relevance to eye disease. Proceedings of the National Acudemy of Sciences, 103(39): 14429-14434 2006.
[40] Gideon Schwarz al. Estimating the dimension of model. The annals of statistics, 6(2): 461-464, 1978.
[41] David Spiegelhalter, Andrew Thomas, Nicky Best, and Wally Gilks. BUGS 0.5: Bayesian inference using gibbs sampling manual (version ii). MRC Biostatisties Unit, Institute of Public Health, Cambridge, UK, pages 1-59, 1996.
[42] Will Stephenson, Zachary Frangella, Madeleine Udell, and Tamara Broderick. Can we globally optimize cross-validatie loss? Quasiconvexity in ridge regression. Advances in Neural Information Processing Systems, 34:24352-24364, 2021.
12

---

[43] William E Strawderman. Proper Bayes minimax estimators of the multivariate normal mean. The Annals of Mathematical Statistics, 42(1):385-388 1971.
[44] Shu Yu Tew, Daniel F Schmidt, and Enes Makalic. Sparse horseshoe estimation via expectation- maximisation. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 123-139. Springer, 2022.
[45] Aad W Van der Vaart. Asymptotic statistics, volume Cambridge university press, 2000.
[46] M Varah, James. On the numerical solution ill-conditioned linear systems with applications ill-posed problems. SIAM Journal on Numerical Analysis, 10(2):257-267, 1973.
[47] Boxiang Wang and Hui Zou. Honest leave-one-our cross-validation for estimating post-tuning generalization error. Stat, 10(1))e413, 2021.
[48] ASAR Yasin, Adnan Karaibrahimoğlu, and GENÇ Aştr. Modified ridge regression parameters: comparative Monte Carlo study. Hacettepe Journal of Mathematics and Statistics, 43(5); 2013.
[49] Yongli Zhang and Yuhong Yang. Cross-validation for selecting model selection procedure. Journal of Econometrics, 187(1):95- 112, 2015.
13

---

A EM Procedure for Multiple Means
We show that the Bayesian EM procedure provides sensible estimates of the regularization parameter even in the setting of the normal multiple means problem with known variance o². In this setting, the LOOCV unable provide any guidance on how to choose À due all the information for each regression parameter being concentrated in a single observation. We use the r parameterisation of the hyperparameter, rather than r², as the resulting estimator has an easy to analyse form.
In the normal multiple means model, we are given (y:)3, ~ [ 1), ie.. y is a p-dimensional normally distributed vector with mean B3 and identity covariance matrix. 'he conditional posterior distribution of [3 is:
B, N ((1 -n),, o 10---]] x))
(17)
where 1/(1 + +2). Under this setring, Strawderman 43 proved that if 2 3, then any estimator of the form
(1-r(#lell))
2 1
y
(18)
where 0 S r (2 [lvll?) S 2 and r(-) is non-decreasing, is minimax, i.e., it dominates least-squares. We will now show that our EM procedure not only yields reasonable estimates in this setting, in contrast to LOOCV, but that these estimates are minimax, and hence dominate least-sqaures.
For the normal means model, can obtain closed form solution for the optimum by solving for the stationary point for which T = Tr. with т ~ C+ (0, 1):
arg min (Faal log p(u13, r) - log p(3(r) - log π [T)]]}
arg min
log 2 +
w
+ log(1 (1 -7)}
2
2T
p+ + V12 Sw 2pw + w
5
2(2 + p)
and with w = -E E[B3] = (1 x3 + (1 - k)p, 8 = ||y|| and T Ime
This yields
(Vpl(n 2)pp 8r 8) 2(x 1)23((n 2)p 4) + (n - 1) 2² np + (x 1)2s)
K
+ 2(2+p)
with solution x = (p 2)/s. Plugging this solution into 17 of ( (17 is the form (18 with
we note that the resulting estimator
2
1 =
[1vll))
p p+
2 1
|
|
p + 2
p-2 2
As we have (111 < when 2 6, the EM ridge estimator is minimax in this setting for 2 6.
B Proof of Theorem 3.1
We prove that for sufficiently large n, continuous injective reparameterization of the negative log joint posterior of (5) & is convex when restricted to 7² > This is sufficient, since unimodality preserved by strictly monotone transformations and continuous injective reparameterizatioos.
Specifically, for the presented hierarchical model, the negative log joint posterior up to an additive constant is
n + + 2
logoz + 20 ly X,3 p 2 2a
1
||
2 log 72
202,2 (a + b) log(1 + +
and reparameterising with Φ = B/0, p 1/a and = 1/T and reorganising terms yields
- (n + + 2)logp (p+2 2a) log x + (a + b) log(1 + +2) + 1 |lpy x 1x
14

---

The first three terms can easily be checked to be convex via second derivative test, for which the convexity of the second term is contingent on the condition that + p/2, condition that holds true in our specific scenario with a = b 1/2. For the last two terms, the combined Hessian is of block form [A, B; BT, CJ with
X X + x² Ip, B 2x and C = || Symmetric matrices
of this form are positive definite if A and its Schur complement
C - --B B = / (x"x/X +
1p)
are positive definite, Clearly, A is positive definite. Moreover, for n > 4/(22m/) we have +1
Φ 6'(X"X/2 + - Φ- p"(VE" /x2 +1)
= E/X +I)=
S vvvVV Øx 1x2/((ng (17r))
=((I ((I- VV 1
IIVVT x13x2/(mmm)
S 1|0112 ((6)2/((2mm))
|
where we used the fact that VVT is the orthogonal projection onto the column space of V. The overall inequality implies the required positivity of the Schur complement.
C Derivation of Equation 9
C.1 Derivation of ESN
Here we show that
Σ E
(()) g000
= tr (Cov [3] + Σ E|3,1
j=1
o"tr(AA,)) + 3|
This is rather straightforward proof. We use the fact that given random variable 2;; the expected squared value of is E [+²]] Var[æ]] E[a]²
C.2
Derivation of FSS
Here we show that
Eg Ily X,31121 (10))820] = Ily X + Cov(g]
(19)
We first provide an important fact on the quadratic forms of random variables in Lemma below: Lemma C.1. Let b be a p-dimensional random vector and be a p-dimensional symmetric matrix. If [b] μ and Var(b) = Σ, then E b" Ab] tr(AE) 4² A4,.
Now, we expand the left-hand side of Equation
Eg Ily x 2
E, (y X.3) X
Eg Ty - 2y x,3 XT
y 1 + +
XE|3]
E
x,3
(20)
The use of lemma C.I allows Equation 20 be rewritten as
Eg ||ly x yTy 2 X B| + X Cov(o])
= - 2| tr X Cov(8)
- x t XA;")
15

---

D Solving for the parameter updates (Derivation of Equation I1
Rather than solving a two-dimensional numerical optimization problem (TO). we show that given fixed 7², we can find closed formed solution for o², and vice versa. To start off, we need to find the solution for o² as a function of 7². First, find the negative logarithm of the joint probability distribution of hierarchy (5):
arg min { [i - log X, - log p(8)7?,0²) - log p(o²) log ([>))]}
(21)
o²
Dropping terms that do not depend on o² yiclds:
arg
min
{Eg [ log M, B, log -8) p(o²))}
x²
arg min
+ p 2
log o²
ESS 2x
ESN
log o²
20
d²
n + p +
ESS 20
ESN
arg min
log 0²
(22)
d²
2
20
Solving the above minimization problem involves differentiating the negative logarithm with respect to o² and solving for o² that set the derivative to zero. This gives us:
0 0a²
{
n1 + 2
ESS
ESN
log 2
2x
20
2 + nl +
ESS
ESN
0
202
2(02222
2(02)272
7'ESS + ESN 6²
(23)
(n + p + 2)72
Next, to obtain the M-step updates for the shrinkage parameter 72, we repeat the same procedure find the negative logarithm of the joint probability distribution and remove terms that do not depend on either a² or 72:
arg min
[Eg | log p(y) loog p(8)72,00 - llog p(o²) -oo π(>))}
{
p-+p++ 2
ESS
ESN
P + log(1 + 72) + log 72
= arg min
log
+
(24)
2
20
20
2
2
Substiting the solution for a² 23 into equation (24). yields a Q-function that depends only on -2.
We eliminate the dependency on o² by finding the optimal の² as function - and substitute into the Q-function of (24):
ESN + +SS
)}}
arg min
(1 p) log 72 + +2 log(1 + r²) + (n +
1 log
[n [m+ 2)7
2
22
(25)
Differentiating 25 with respect to 72 and solving for the 7² that set the derivative to yields:
ESN + 7SSS
p) log 72 + log(1 + 2 (n + p log
0r2
(n ++ 2)
(EESS F:SS)) (FSN ESNN + ESS FSSSp)
ESN F.SNn
0.
2- +
(26)
The 7 update the positive solution the quadratic equation (in terms of 72) 26
#2 (n - 1JESN - (1 + PSS + /(4n + 4)ESN(3 + ESS + ((1 n) + (p + 1)SSS)
(6 ESS
16

---

Table 3: Pseudo-code of EM algorithm with complexity of individual steps.
EM Algorithm with SVD
Operations
Input:
Standardised predictors X E R=*x, centered targets y E IR" and convergence threshold 6
Output: B E TRP
pi IF 2 N
O(1)
[U. Σ, V sx)(X)
O(mr))
s² [ 1,...., D2
O(r)
c= (U O)@s
os
O(nr)
ELSE
[V. Σ eigen X X
O (mr?)
s² = [
((r)
== V'x
(np)
y-y'y
(
72
0(1)
o² (I/n) Ei=1 (1i (11-0)² 9 -1/n)) Ei 1i
O(n)
RSS + 00
0(1)
DO
RSS t RSS
O(k)
0j $4 1/72
Cj
0(kr)
(E-step)
ESN + Σ 2 + t
Σ
+ max(p - n,
0(kr)
+ +
j=1
RSS +- Y -2 Σ
0jCj
O(kr)
5=1
j=1
ESS < RSS 0²
O(kr)
+
(M-step)
+ (4n + 4)EESN (3 + EJSSS + ((1 n)ESN + + 1)ESS)?
O(k)
72 (n 1)EN (1 F)SS v/g
O(k)
(6 + 2p)
0² EESS ESN
O(k)
(n + p + 2)
RSS RSS|
O(k)
(1+ (RSS)
until 8 <E
02 87 + 1/72
Cj
O(kr)
3= Va
O(pr)
return |3
17

---

Table 4: Pseudocode of the fast LOOCV algorithm with complexity of individual steps. R has column vectors r for 1 S j r.
Fast LOOCV ridge with SVD
Operation
Input:
Standardised predictors X E Ruxp, centered targets E IR" and a grid of penalty parameters I. ..)
Output: / E IRP
pi IF p 2 n
O(1)
[U, Σ, VJ = svd(X) (E1,1,..., Enr)
O
O(r)
R (stU)...., S,
((r)
yU")@ Os
O(nr)
ELSE
[V. Σ eigen X
0
s² E 1'
E
O(r)
R= XV
(nrp)
c R'y
((r))
U (r1/s1, kr//r)
for X E L
h4 Σ
j
unz.
(i= 1.... n.)
8) +
j=1
0j Cj
0
} +
e= y Ra
0
Ei
CVE(A)
O(In)
Find X" arg min (CVE(A)) det
0(1)
Cj
0j 5 X
0(r)
G Va
O(pr)
rettur
18

---

E Supplementary Results Material
E.1 Real Datasets Details
Table 5: Real datasets details
DATASETS
ABBREVIATION
71
P
TARGET VARIABLE
SOURCE
BUZZ IN SOCIAL MEDIA (TWITTER)
TWITTER
583250
77
mean number active discus-
UCI
sion
BLOG FEEDBACK
BLOG
60021
281
number of comments in the next
UCI
24 hours
RELATIVE LOCATION OF SLICES
CT SLICES
53500
386
reference: Relative image loca-
UCI
ON AXIAL AXIS
tion on axial axis
BUZZ IN SOCIAL MEDIA (TOM'S
TOMSHW
28179
97
Mean Number of display
UCI
HARDWARE)
CONDITION-BASER MAINTENANCE
NPD COM
11934
16
GT Compressor decay state efficient
UCI
NAVAL PROPULSION PLANTS
CONDITION-BASED MAINTENANCE OF
NPD TUR
11934
16
GT Turbine decay state coeffi- cient
UCI
NAVAL PROPULSION PLANTS
PARKINSON'S G
PT MOTOR
5875
26
motor UPDRS score
UCI
PARKINSON'S TELEMONITORING
TOTAL
5875
total UPDRS score Rings (age in years)
UCI
ABALONE
ABALONE
4177
UCI
COMMUNITIES AND CRIMB
CRIME
1994
128
ViolenCrrimeserPop
UCI
AIRFOIL SELF-NOISE
AIRFOIL
1503
Scaled sound pressure level (decibels)
UCI
STUDENT PERFORMANCE
STUDENT
649
grade (with G2 re- moved)
UCI
CONCRETE COMPRESSIVE STRENGTH
CONCRETE
1030
Concrete compressive strength (MPa)
UCI
FOREST FIRES
F.FIRES
517
forest burned area (in ha) Median value of owner- occupied homes in $1000's Total Interactions (with com- ment, like, and share columns removed)
UCI
BOSTON HOUSING
HOUSSIGG
506
13
FACEROOK METRICS
FACEBOOK
500
19
UCI
DIABETES
DIABETES
442
10
quantitative measure of disease progression one after base- line
13
REAL ESTATE VALUATION
E.ESTATE A.MPG
414
house price of unit area city-cycle fuel consumption in miles per gallon
UCI
AUTO MPG
398
UCI
YACHT HYDRODYNAMICS
YACHT
308
7
residuary resistance per unit weight of displacement price
UCI
AUTOMOBILE
A.MOBILE
205
UCI
RAT EYE TISSUES
EYE
120
200
the expression level of TRIM32 gene
RIBOFLAVIN
RIRO
71
4088
log-transformed riboflavin pro- duction
CROP
CROP
24000
3072
24 crop classes
UCR
ELECTRIC DEVICES
ELECD
16637
4096
electric devices
UCR
STARLIGHT CURVES
STARL
9236
7168
starlight curves
UCR
19

---

1.0
0.8
CLM Grid
0.6
0.4
0.2
1.0
CV Fixed Grid
0.0
0.4
0.5
1.0
0.4
0.6
1.0
BayesEM
BayesEM
BayesEM
Figure 3: Comparison of predictive performance (R²) of EM algorithm (r-axes) against CV with fixed grid (y-axes, top) and glmnet heuristic (y-axis, bottom). Columns correspond to the results of linear features (left), second-order features (middle), and third-order features (right). Negative values are capped Points skewing toward the bottom right indicate when our EM approach is giving better/same prediction performance as LOOCV (colored in green).
20


Reviews:
[{"review_id": "6T6SeP716P", "reviewer": "Reviewer_Ahob", "summary": "This paper proposes a new method for estimating the regularization parameter in ridge regression models that is faster, than the widely used leave one out cross validation (LOOCV) technique. The new method is derived from a Bayesian view for estimating this regularization parameter which, through standard Expectation Maximization procedure, lends itself to a nice optimization problem with two variables. The efficiency comes from the approach taken to solve this optimization problem; the Authors are able to formulate the solution using exact closed form expressions, together with a pre-processing step -- SVD for the data matrix -- that is computed only once at the beginning of the procedure. The paper demonstrates the accuracy and efficiency of the proposed method using an extensive set of experiments on simulated\nand real world datasets, and compares its performance against the fast LOOCV with a fixed grid for the ridge/regularization parameter, as well as with the grid based on data-dependent heuristic available in R and Python packages for glmnet.", "questions": "No questions at this point.", "limitations": "See my remark on missing discussions and point that require further elaboration.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The contribution of this paper is clear and interesting, and the derivation for the procedure is straight forward to follow. The paper is nicely written, well organized, and the experiments are extensive and sufficient. I believe this will be a good addition to the literature.", "weaknesses": "The missing pieces that require discussion and elaboration from the Authors is when this proposed method will be worse than the fast LOOCV? Any assumptions on the data distribution? And how sensitive is this method to the choice of $a$ and $b$ in Eq. 9?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback and especially your assessment that the paper will be a good addition to the literature. Regarding the questions you identified:\n\n**When is the proposed method worse than the fast LOOCV?**\n\nWhile this is a natural question, it is from our perspective not easy to answer. In contrast to the LOOCV estimator, which has an obvious structural weakness in sparse covariate settings where every data point carries important information, the EM does not exhibit such an apparent weakness. This may be seen in the setting in which a covariate, say $x_i$ is non-zero for only a small number of data points, say $n_1 \\ll n$ -- a situation common in problems such as genomics. In this case, the information contained in the data about coefficient $\\beta_i$ is concentrated in this small number of data points -- the remaining data $n-n_1$ datapoints offer no information about $\\beta_i$. In this case, LOOCV will exhibit a high variance as dropping one of the $n_1$ relevant datapoints will have a substantial impact on the resulting estimate of $\\beta_i$. In the extreme case that there is a single non-zero value per covariate, i.e., $X=I_n$ (the well studied multiple means problem, which is representative of a range of real world problems including orthogonal function smoothing) the LOOCV procedure offers no way of estimating $\\lambda$ as each datapoint we drop removes *all* the information about one of the coefficients.\n\nIn contrast, the Ridge EM approach, by exploiting the Bayesian interpretation of the coefficients as being generated by an exchangable random (Gaussian) process is not only able to provide estimates of $\\lambda$ in this setting, but is actually shown to yield estimates that are minimax, and therefore strictly superior to least-squares, as shown in the Appendix.\n\n**Any assumptions on the data distribution?**\n\nThere are no sensitive assumptions on the data distribution given that we employ a weakly informative prior. In particular, for large $n$ the EM algorithm is guaranteed to target a simple unimodal posterior distribution (see [global rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq&noteId=OV1Ci3vVUB)). In this regime, it requires few iterations to converge to a minimum resulting in the good documented time complexity. While we are lacking a theoretical result for small $n$, we empirically observe in this regime, in a wide range of settings, an advantage of EM over LOOCV in terms of prediction risk.\n\n**How sensitive is this method to the choice of $a$ and $b$**\n\nThe choice of $a$ and $b$ does not significantly impact our method's performance because the class of priors used in the paper has a heavy tail. Specifically, we select a half-Cauchy prior ($a=1/2$, $b=1/2$) as it is well-studied and recommended as a favorable default choice for a global shrinkage parameter. See also the [global rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq&noteId=OV1Ci3vVUB).\n"}, {"review_id": "Y4ib4qyPne", "reviewer": "Reviewer_nvLc", "summary": "This paper studies the comparison of cross validation (CV) and  expectation maximization (EM) on tuning the regularization hyper-parameter $\\lambda$.\nThe authors explained the running time of these two methods and showed that EM is faster than CV.\nThen, the authors also provided some experimental results to support the conclusion that EM is faster than CV and demonstrate that they have comparable predictive performances.", "questions": "See weaknesses", "limitations": "na", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The proposed EM approach for tuning $\\lambda$ is quite interesting.\nI have not seen this approach before.", "weaknesses": "The result is mostly experimental.\nMost of the theoretical analyses are on the running time.\nThe authors may want to provide some theoretical explanation on why these two methods end up having comparable predictive performances.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive evaluation. Regarding some further theoretical insights we refer to the [global rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq&noteId=OV1Ci3vVUB). \n\nWhile we do not currently have a specific explanation for the typically similar predictive performance, we have examined more closely the cases in which Ridge-EM is expected to be superior to LOOCV; see [rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq&noteId=es7IWesb9i)."}, {"review_id": "Jf1RJKCgIk", "reviewer": "Reviewer_iDip", "summary": "This article proposes a new method based on an Expectation-Maximisation framework to train a Bayesian formulation of ridge regression for penalised regression. It is compared with the standard cross-validation approach for training the regularisation parameter lambda, instead iteratively updating the value algorithmically. Novel data pre-processing via SVD is proposed that accelerates not only the EM algorithm but also the standard glmnet algorithm. Extensive comparisons are made between the proposed EM method and the now-popular CV approach and computational costs are considered.", "questions": "What is stopping this from becoming the norm if it is superior to glmnet/scikit-learn? One thing that people like about glmnet (not necessarily a good thing in statistical terms) is that it spits out specific predictors with specific values that are easy to interpret. Does this offer an end product that is similarly easy to interpret without any real domain expertise in inference to get it work? If so, is it necessarily a good thing to let non-statisticians use this method without understanding what's happening under the hood?\n\nWhat is the behaviour in contrast to CV with potentially non-i.i.d. Data? Does the new method make an analogous assumption to LOOCV, or avoid it? Can this difference be identified empirically?\n\ncv.glmnet solves for lambda.min and lambda.1se, i.e. the lambda value that minimises the objective function vs. one that is inflated a little for a little extra sparsity. Which is being compared with in the results here? And is there an analogous decision to make with the EM method? \n\nAre people likely to object to using explicit priors in a Bayesian method vs regularisation for a frequentist method? Are you able to convince people of the use of why it’s acceptable to use the ones you have chosen generically for any data set? What is the sensitivity to pi(tau)?\n\nIf you were really doing things at scale and worried about computational costs, it would be fairly straightforward to distribute the fixed grid evaluations of LOOCV in parallel with minimal overhead, effectively removing the factor of L that disadvantages the method. To what extent will the Bayes-EM formulation take advantage of similar distributed computing if available, and to what extent does this undermine the performance advantage represented here?\n\nAre there specific instances in which the EM algorithm doesn’t converge quickly? Are there good diagnostics for identifying when this is happening? Figure 2 seems to show slightly concerning behaviour in terms of number of k needed for increasing p. Do you have some theoretical characterisation of this bump, possibly representing the degree of non-Gaussianity of the posterior? It would be nice to know a priori what to expect from the EM algorithm, so a slow convergence rate will not come as a surprise. Presumably in this domain the LOOCV will have the advantage in a performance-time trade-off here, which is not super clear from the results as presented.\n\nLOOCV basically picks more or less at random from a pool of highly correlated predictors, which is not a good thing. Does BayesEM exhibit the same behaviour? Does this Bayesian formulation open up the possibility of adequately modelling correlation of coefficients in sparse regression via modelling some correlated beta posteriors?\n\nHow robust is this to non-normally distributed betas? One advantage of LOOCV is that the regularisation is an inferential decision rather than a modelling decision, in contrast to BayesEM, where you are making positive belief statements, which may be empirically false. If they are false, how much does that actually matter?\n", "limitations": "No immediate negative societal impact to be concerned about. It seems to be a reformulation of a method that is already in widespread use.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I enjoyed reading this paper. The very widely used model specification of penalised regression is reformulated in a Bayesian framework and the algorithm approaches the problem in a principled manner. The advantages compared to the standard cv.glmnet or python equivalent are fairly convincing. The algorithmic steps are explained clearly and the results are impressively thorough (24 data sets!). The usability of the method appears to be high, with the practical hurdles for using this in comparison with glmnet or similar appear to be low: this seems like it could become competitive with the present standard for biomarker analysis if presented and promoted appropriately, although the same could also be said of many other methods that are now in widespread use presently.", "weaknesses": "Some of the phrasing of the results is quite unclear and possibly suspicious: “Up to one order of magnitude” - that can mean a factor of two or a factor of ten, depending on who you ask.  Also “up to” suggests that this is the maximum difference rather than a bulk estimate that might be more informative in general. Similarly, towards the end: “Up to 49 times faster” - I would be more convinced by a summary involving a characterisation of average differences rather than a maximum value, that might be an outlier or a particularly flattering data set. If that average needs to be characterised as a function of p or n or similar to be accurate then so much the better.\n\nSome result cells are empty in Table 2: it is OK for these to not be present, but it is not explained why, so explaining more clearly to avoid suspicion would be nice.\n\nIt seems that some of the inspiration for this method comes from relatively recent related work: “The EM procedure discussed in this paper is based on a recently published procedure for learning sparse linear regression models” - please be more clear about the specific contribution of this work vs previous work, i.e. horseshoe vs ridge or similar.\n\nSome of the language and explanation could be tidied up slightly. Some examples:\n\n- \"that minimizes (3.” - typo?\n- There is no real difference between “missing data” and a latent variable in a Bayesian framework. If you could explain more clearly what insight is meant to be provided then that would be nice: if there is a clever trick here then please spell it out for those of us who haven't seen it before.\n- “Described in 5”? “Described in Section 5” maybe?\n- I assume the third order models include linear/quadratic terms and similar, but it’s slightly unclear.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewers' positive comments on our method's novelty, high usability, and principled presentation, and we are delighted to learn that the reviewer enjoyed reading the paper.\n\nImportantly, we want to emphasize that both the EM method and LOOCV share the same underlying model, which is ridge regression. This model penalizes the sums of squares and assumes a normal distribution over the regression coefficients (beta). The assumption of normality over beta is explicit in the Bayesian treatment of ridge regression, where explicit priors are employed, and it is implicit in any form of ridge regression, making the underlying model targeted by both methods identical. The primary distinction between the two methods lies in how they estimate the optimal lambda. With this context in mind, we address specific questions:\n\n**How robust is this to non-normally distributed betas?**\n\nThe robustness of a method to non-normally distributed betas depends on the underlying model chosen, such as ridge regression in this case. Given that ridge regression is not necessarily the best choice when the betas deviate from a normal distribution, neither the EM method nor LOOCV will be optimal in this situation. However, we note that ridge regression may still be calibrated to dominate unpenalized least-squares for regardless of the distribution of the coefficients as long as the resulting estimator is minimax. Our results regarding minimaxity suggest that our Ridge EM procedure is calibrated in this fashion.\n\n**LOOCV picking random variables from a pool of highly correlated predictors …..  Does this Bayesian formulation open up the possibility of …. sparse regression…….?**\n\nWe would like to clarify that this issue is more commonly associated with Lasso regression rather than ridge regression. Ridge regression does not perform variable selection or induce sparsity, while Lasso does. To address variable selection with highly correlated predictors, one may opt for elastic net regression. LOOCV can be used to find the optimal lambda parameter in elastic net regression as well. Similarly, our EM algorithm can also be applied for optimizing lambda in elastic net regression, but further studies are required to derive efficient updates and conduct a thorough comparison of computational efficiency and prediction performance between the two methods. So the answer is yes, there are Bayesian approaches that adequately address this issue but it is important to note that our current paper focuses solely on ridge regression, and the exploration of Bayesian methods for variable selection in this context lies beyond the scope of this paper.\n\n**What is the behavior in contrast to CV with potentially non-i.i.d. Data? Does the new method make an analogous assumption to LOOCV, or avoid it?**\n\nFrom a high-level perspective, we expect our method's sensitivity to non-independent data to be similar to LOOCV, as both target the same underlying model. However, when the correlation structure is known, the EM procedure can easily be adapted.The process itself remains minimally changed; specifically, we perform an SVD on L*X, where L represents the Cholesky decomposition of the covariance matrix, while retaining the rest of the methodology intact. In contrast we are uncertain about the adaptation process of LOOCV to fit the context of correlated data. While the possibility of a straightforward adaptation exists, we currently unaware of this.\n\n**cv.glmnet solves for lambda.min and lambda.1se…... Which is being compared…..?**\n\nWe opt to compare against lambda.min as selecting lambda.1se seems inappropriate for ridge regression, given its inability to induce sparsity. Currently, our answer to the latter question is no; however, a possible approach might involve considering the posterior variance to assess a similar version of lambda.1se. Nonetheless, efficiently evaluating this approach remains an open problem and warrants further investigation.\n\n**Use of multiple cores to speed up CV**\n\nIn the context of a finite $n$ situation, we concur with the reviewer's observation that distributing fixed grid evaluations of LOOCV in parallel with minimal overhead is a valid and valuable approach. However, it is essential to consider that in practical scenarios or when dealing with increasing sample sizes (n), the overhead for parallel processing might become substantial. As $n$ increases asymptotically, the need to work with larger grid sizes may eventually exceed the capabilities of any finite computing architecture.\n\n**Some result cells are empty in Table 2……. explaining more clearly……..**\n\nThank you for bringing this to our attention; we will ensure to provide a clear explanation in the text. Our original plan was to include all results in the table. However, due to time constraints, we had to set a maximum run-time of 3 hours for each dataset, and any settings that exceeded this time limit were subsequently excluded from the table.\n\n**I assume the third-order models include linear/quadratic terms and similar, but it’s slightly unclear.**\n\nThe third-order models include 3-way interactions, 2-way interactions, squares, and cubics of the predictors. \n\n**What is the sensitivity to pi(tau)?**\n\nThe sensitivity is very low due to the robustness of the chosen prior (see Global rebuttal for more details).\n\n**please be more clear about the specific contribution of this work vs previous work, i.e. horseshoe vs ridge or similar.**\n\nIn this work, we show how the specific structure of ridge regression can be exploited to derive substantially faster EM updates in this important special case (the updates for the general case are O(p^3) vs O(p) that we show is possible for ridge regression). Additionally, we have demonstrated how to use similar techniques to speed-up LOOCV and we have provided a comprehensive comparison of both approaches in terms of computation time and predictive performance. None of these questions have been addressed by prior work."}, {"review_id": "WcSOHqWi8s", "reviewer": "Reviewer_3fTV", "summary": "This paper proposes a simplified implementation of the EM algorithm of the Bayesian Ridge regression for the hyperparameter estimation. The advantage of the proposed method is in its computational cost when searching the hyperparameter. The numerical experiments support their proposal in terms of computational time. \nThe essential idea is to use SVD and compute all necessary quantities at preprocessing, and to lighten the computational load of hyperparameter search in the main loop by using them.\n", "questions": "- What $k$ and $q$? Please explain about them when they first appear in the text. \n", "limitations": "The authors address their research limitations well. I think there is no concern about the potential societal impact. \n", "rating": 5, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "Originality & Significance: This kind of trick (completing as many of the necessary calculations as possible as preprocessing) is, in my opinion, done by many researchers actually, though it may not be mentioned much in papers or textbooks. Especially, using SVD in preprocessing is a well-known technique in the family of linear models. It may be new to point out that using EM via Bayesian is less computationally expensive than LOOCV with standard ridge regression in the order of $n$, but I do not think this is a special property of ridge regression and is not expected in more general settings. For these reasons, I consider the originality and significance to be low. \n\nQuality: All the numerical experiments seem to be certainly conducted.　Numerical codes are also distributed.　In this sense, the quality is high.\n\nClarity: Most of the contents are well written, but some important parameters characterizing the computational costs such as $k$ and $q$ are not well explained. For this reason, the clarity is at the normal level.\n", "weaknesses": "- The basic idea of the research is not surprising. \n\n- The Bayesian ridge regression is different from the ridge regression in the level of concept, although the estimators tend to behave similarly. In that sense, I am not sure if it is suitable to compare these two things (EM in Bayes and LOOCV in the standard ridge.). \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the positive comments on the simplicity and effectiveness of the proposed method. However, we disagree with the assessment of our work’s significance. In particular, our technical contributions extend beyond merely proposing a singular-value decomposition (SVD) as pre-processing. \n\nUsing an SVD serves only as the starting point based on which we address the challenging questions surrounding the efficient computation of the conditional expectations of regression coefficients and the sum of squared errors. While Sec. 4.1 outlines the necessary steps in a succinct way, this consolidated presentation should not suggest that they are straightforward to derive.  Indeed, the arguably simpler derivations that we also provide for the case of fast LOOCV appear not commonly known (they are not implemented in `glmnet` and `scikit-learn`) despite the fact that this estimation method has been widely used for decades.\n\nRegarding the reviewer's concern about comparing LOOCV to the EM method, we note that, by focussing on the posterior mean of the regression parameters, there is no practical difference between the classical and the Bayesian description of ridge regression—the latter only describes the underlying assumptions explicitly in the form of a normal prior. What we are left with are two different ways to estimate the optimal $\\lambda$, and comparing these in terms of the prediction risk and computation times addresses the valid and important question, which one should be preferred in practice.\n\nResponse to the question:\n\n**What are $k$ and $q$? Please explain about them when they first appear in the text.**\n\nWe are happy to revise the text to have these definitions stand out more: $k$ denotes the number of iterations required for the proposed EM algorithm to converge, and $q$ indicates the number of target variables in a multitarget regression model.\n"}, {"review_id": "ZDiFO70yeX", "reviewer": "Reviewer_fBDz", "summary": "This paper tackles the problem of choosing the regularization parameter in ridge regression. The authors introduce a Bayesian model that has the ridge estimate as its posterior mode with a prior distribution over the regularization parameter; they then propose an EM algorithm to estimate the optimal regularization parameter. The authors compare their estimator to leave-one-out CV (LOOCV), and note that each iteration of their EM algorithm is faster than evaluating LOOCV for one value of the regularization parameter. In experiments, the number of EM iterations is modest, so the author's algorithm is faster than LOOCV.", "questions": "Overall, I think I'm feeling like the current paper has enough novelty in it for publication. But if the authors could argue against this, that'd be the most likely way to change my score -- I'm also curious what the other reviewers think on this issue.", "limitations": "I don't think limitations of the proposed method have been addressed in that they haven't really been studied -- it would be surprising for the proposed EM method to have no failure modes at all, but none are studied in the paper. \n\nI don't think there is any potential negative societal impact of this work.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This paper is well-written and easy to follow (maybe in the top 50% of NeurIPS papers in this regard). I thought the authors made their contributions very clear, and they have extensive experiments. I also think the authors did a good job of backing up all of their specific claims with either experiments or theory.", "weaknesses": "**Editing this here on 8-24-2023** because it seems I can't list authors as readers to any new replies at this point.\n\nThe global comment proving the unimodality of the EM objective is great! I'm convinced, and will raise my score. The strength that I see of this paper is that it is concretely improving on LOOCV in two ways: it has a runtime improvement per iteration (and solid empirical evidence that the number of iterations required for convergence in the EM algorithm is small), and it is guaranteed to be unimodal as soon as $n$ is modest compared to the smallest eigenvalue of the covariance matrix of the covariates. A couple remaining notes that I think could be good to address in a potential camera-ready:\n\n1. LOOCV can be guaranteed to have only one local minimum -- see Stephenson et al. (2021). But, compared to the conditions under which the current authors' EM objective is guaranteed to have one minimum, the conditions in Stephenson et al. (2021) are way more complicated, and their proof is way more complicated. And their experiments show that LOOCV will often have multiple local optima in normal conditions. I think these results make the authors' result in the global comment even stronger.\n2. I still think the paper could be made stronger by an empirical comparison to gradient descent on the LOOCV objective. But I think the paper is solid as-is with the inclusion of all the discussion above and the new proof in the global comment.\nBest, fBDz\n\nWilliam T. Stephenson, Zachary Frangella, Madeleine Udell, Tamara Broderick. Can we globally optimize cross-validation loss? Quasiconvexity in ridge regression. NeurIPS 2021.\n\n---------------------------------------------------------------\n\n\nI think a method that in both accuracy and computation beats LOOCV for ridge regression would be big news for the ML community. I'm not entirely convinced that the proposed method is doing that, though. And something that is big news also needs to come with a lot backing it up. First, a lot of the issues with LOOCV in the paper seem to stem from the choice of a finite grid over the regularization parameter -- the authors point out that it is slow to evaluate LOOCV for all points in the grid, and also the choice of grid can be tricky. It seems like the authors' method only gets around these issues by virtue of being an optimization algorithm. However, there is a much more common optimization algorithm: gradient descent over the LOOCV objective. True, the authors point out that the LOOCV objective can have local minima. But it often doesn't. How do we know the same isn't true of the authors' method? I really think a detailed comparison to this approach is needed, as gradient descent seems to have many of the same strengths as the authors' EM method and an equal number of weaknesses (i.e. I think we should assume both can have issues with local optima until proven otherwise).\n\nMy second main issue with the paper is that it doesn't really justify where this EM algorithm came from or what it's really doing. **Why** is this algorithm outperforming LOOCV? What is it picking up on that LOOCV is not? Why not use some other miscellaneous approach? E.g. could we sample from the posterior and take the posterior mean of $\\tau$? How do the chosen priors on $\\sigma$ and $\\tau$ affect things, and why are those good priors on $\\sigma$ and $\\tau$? Overall I this approach needs to be a lot more explained and justified if the proposal is to replace LOOCV with it.\n\nOverall, I think the current contributions of the paper aren't quite enough. The derivation of the fast LOOCV algorithm is straightforward linear algebra (I agree it's surprising that scikit-learn isn't implementing this method, but I think that's just more of a criticism of scikit learn). And the derivation of the EM algorithm iterates seems relatively straightforward from [26]. I think the real-data experiments are pretty good, but I'm not sure this is enough to carry the paper. It would be really interesting to investigate the proposed EM algorithm a bit more -- e.g. get into the motivation as above or study it theoretically (e.g. under what conditions is the assumption on line 208 that the number of iterations is $o(n)$ correct?)\n\nI had a few other more contained comments:\n- I don't think $q$ was defined mathematically in the paper.\n- On line 90, some non-standard notation for inner-products is used that I don't think was defined (maybe $\\tilde x_i^T \\hat\\beta_{(i)}$ was meant?)\n- The authors criticize LOOCV by saying \"For the squared-error risk to be minimax, and therefore donimate least-squares ... the number of grid points [in LOOCV] $L$ cannot be finite ...\". Does the proposed EM algorithm fix this problem? In particular, does it converge in a reasonable number of finite iterations, and, if not, is it still minimax if it is not run for infinite iterations?\n- In equation (8), what the expectations and covariances are with respect to should be specified -- I don't think any data-generating distribution was specified in the paper.\n- Line 147 \"we provide the derivations of (8) in the Appendix\" -- it should be specified where exactly in the appendix this can be found. This comment is true of every reference to the appendix throughout the paper -- readers shouldn't have to hunt through the whole thing to learn more about a specific point!\n- Line 174: is $s_j$ the same as $\\Sigma_{jj}$? If so, notation should be standardized. If not, I didn't see where $s_j$ was defined.\n-In equation (12), it should be clarified that the $\\tau^2 (p - n)_+$ is outside the summation.\n- Figure 2 is about showing that LOOCV is much slower than the EM algorithm. But the actual difference in runtime is like half a second at its maximum. I think it would be much more compelling to show a difference in runtime that is computationally meaningful.\n- Line 218: \"find the smallest value of $\\lambda_{max}$ that sets all the regression coefficient vector to zero\" -- isn't this typically $\\lambda_{max} = \\infty$ for ridge regression? \n\nI have a few other", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for the positive comments and the assessment that the paper is likely impactful. \n\n**Regarding your concern about novelty:** We show for the first time how the specific structure of ridge regression can be exploited to derive substantially faster EM updates in this important special case. Additionally, we show how to use similar techniques to speed-up LOOCV and we have provided a comprehensive comparison of both approaches in terms of computation time and predictive performance. None of these questions have been addressed by prior work. Technically, one should also not underestimate the derivation of the EM algorithm iterates. Even though Sec 4.1 presents a concise outline of the necessary steps, their derivation is far from straightforward.  Indeed, the arguably simpler derivations that we also provide for the case of fast LOOCV appear to be not commonly known (they are not implemented in `glmnet` and `scikit-learn`) despite the fact that this estimation method has been widely used for decades.\n\n**Regarding the idea of using gradient descent as an alternative to avoid pre-specified grid for LOOCV**: We agree that this is an interesting direction to explore. However, note that this would require finding an appropriate schedule of the step size, which might be non-trivial to do and, if chosen poorly, might cause the approach to be inferior to a fixed grid. Moreover, each evaluation of the score and the gradient still incurs $O(mn)$ operations, resulting in computational overhead compared to the EM procedure, which requires only $O(m)$ operations per iteration, where $m$ is ${\\rm min}(n,p)$.\n\nIn terms of the convexity of the optimization problem, we demonstrate that as $n$ approaches infinity, the posterior distribution becomes unimodal thus ensuring convexity (see [global rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq&noteId=OV1Ci3vVUB)). However, we cannot guarantee the same for the LOOCV.\n\nResponse to specific questions:\n\n**Why is this algorithm outperforming LOOCV?**\n\nThe EM approach offers a key advantage over LOOCV since it avoids omitting any data point during model fitting, making it well-suited for sparse covariate settings commonly seen in genomics, where coefficient information is concentrated in just a few observations. In extreme cases like the multiple means problem (${\\bf X} = I_n$) [15], where each coefficient’s information concentrates in a single observation, LOOCV becomes uninformative for $\\lambda$ selection, while RidgeEM explicitly connects coefficients through a probabilistic Bayesian interpretation of $\\lambda$ as the inverse-variance of the unknown coefficient vector (more details in lines 108-119 of the submission). In non-sparse covariate settings, both EM and LOOCV demonstrate similar prediction accuracy (Section 6.2), but we outperform in terms of computational complexity due to the reduced main loop complexity as summarized in Table 1.\n\n**Why not use some other miscellaneous approach? (E.g. could we sample from the posterior and take the posterior mean of $\\tau$)?**\n \nIndeed, we can estimate using the posterior mean obtained through sampling from the posterior distribution. The primary concern with this is efficiency, as sampling introduces higher time complexity compared to EM and LOOCV. \n\n**How do the chosen priors on $\\tau$ and $\\sigma$ affect things, and why are those good priors?**\n\nWe note that these priors have been chosen because they are standard choices that are weakly informative/uninformative and result in a high degree of posterior robustness (see the global rebuttal for more details).\n\n**Line 218: \"find the smallest value of $\\lambda_{max}$ that sets all the regression coefficient vector to zero\" -- isn't this typically $\\lambda_{max} = \\infty$  for ridge regression?**\n\nThank you for pointing this out. We will correct this in the paper. We adopt the approach used in the glmnet package in R, selecting the largest  value for when $\\alpha=0.001$. The penalty function used by glmnet is $\\lambda[(1-\\alpha)||\\beta||_2^2 + \\alpha||\\beta||_1]$ with $\\alpha=0$ corresponding to ridge regression. The sequence of $\\lambda$ is derived from this with  $\\alpha=0.001$. This underscores the challenges in specifying a grid for LOOCV methods, as we believe there are potentially several other approaches to handle it, lacking consistency and sufficient evidence to support their ideas.\n\n**In equation (8), what the expectations and covariances are with respect to…..?.**\n\nThe conditional expectations and covariances are computed with respect to $\\beta$. These follow a multivariate normal distribution as specified in hierarchy (4)-(5).\n\n**The authors criticize LOOCV by saying \"For the squared-error risk to be minimax, and therefore dominate least-squares ... the number of grid points [in LOOCV] L cannot be finite ...\". Does the proposed EM algorithm fix this problem? In particular, does it converge in a reasonable number of finite iterations, and, if not, is it still minimax if it is not run for infinite iterations?**\n\nWe acknowledge that in a finite number of iterations, we may not attain exact minimaxity. However, we note that we can attain some form of $\\epsilon$-minimaxity by allowing the algorithm to run for longer. In contrast, even if we allow the grid size to be arbitrarily fine (a necessary condition), we are not aware of any proof that the minimaxity of the LOOCV estimate will actually be minimax.\n\n**Thank you for bringing the notational inconsistencies to our attention. Below are the responses addressing all questions related to the notational issue:**\n- $q$ indicates the number of target variables in a multitarget regression model\n- On line 90, we intended to express $\\hat{y_i} = x_i^T \\hat{\\beta}(\\lambda)$, where $\\hat{\\beta}(\\lambda)$ represents the estimated regression coefficient when using a specific $\\lambda$ value to predict $y_i$.\n- $s_j =\\Sigma_{jj}$, with $s^2_j$ being the eigenvalue for $X^T X$ as described in line 174."}]

Author rebuttal:
We thank all reviewers for their thoughtful comments and mostly positive evaluation of our work. Some reviewers encouraged us to provide further theoretical insights into the behaviour of the proposed EM method. Here, we summarise some new considerations that we intend to include in the final version.

**Convergence/unimodality of posterior**

Several reviewers enquired regarding convergence properties and the solution found by the EM algorithm. We note that EM algorithms are generally guaranteed to converge to an exact mode of the posterior, but the solution may be a local mode. While we do not currently have a proof of unimodality for finite sample sizes, we note that the posterior density converges to normality as the sample size tends to infinity. This can be established using the Bernstein-Von Mises theorem (see [1], Section 10.2), as our problem satisfies the relevant conditions: both the Gaussian-linear model $p(y | \beta,\sigma^2)$ and the marginal distribution $\int p(y | \beta,\sigma^2) \pi(\beta | \tau^2) d\beta = p(y | \tau^2,\sigma^2)$ are identifiable, have well defined Fisher information matrices, and the priors over $\beta$ and $\tau$ are absolutely continuous.

**Choice of priors**

Several reviewers asked regarding the sensitivity of the procedure to the choice of the prior over the hyperparameter $\tau^2$, and the specific choice of hyperparameters for the prior ($a=1/2, b=1/2$, yielding a half-Cauchy). We first note that the half-Cauchy is a heavy tailed, weakly informative prior that is frequently recommended as a default choice for scale-type parameters (such as $\tau$), and has been shown to have good properties when used as a prior for ridge-type hyperparameters (see [2]). 

Further, we note that it is very insensitive to the choice of $a$ or $b$. From Theorem 6.1 in [3] we see that the marginal prior density over $\beta$, $\int \pi(\beta | \tau^2) \pi(\tau^2 | a,b) d\tau^2 = \pi(\beta | a,b)$ has polynomial tails in $||\beta||^2$ for all $a>0, b>0$, and has Cauchy or heavier tails for $b<1/2$. This type of polynomial tailed prior distribution over the norm of the coefficients is highly robust/insensitive to the overall scale of the coefficients, which is likely to be *a priori* unknown. We note that this is in contrast to other choices of prior distributions for $\tau^2$ such as the inverse-gamma which are often used but lack this robustness property and are highly sensitive to the choice of hyperparameters.

The prior for $\sigma^2$ is the standard scale-invariant prior that is uninformative and expresses *a priori* ignorance regarding the scale of the data. This is a standard choice.

[1] A v.d. Vaart, "Asympotic Statistics".

[2] N.G.Polson, J.G.Scott, “On the Half-Cauchy Prior for a Global Scale Parameter”, Bayesian Analysis, 2012

[3] O. Barndorff-Nielsen, J. Kent, M. Sorensen, “Normal Variance-Mean Mixtures and z-Distributions”, International Statistical Review, 1982


Meta-review: This paper presents a nice contribution towards the literature on tuning regularization parameters in ridge regression. While ridge regression is often considered a very simple model, the presented results are crisp and clarify some fundamental questions. Three of the reviewers argued rather strongly in the discussion for the acceptance and the area chair is of the same opinion.
Decision: Accept (poster)