Bayes beats Cross Validation: Efficient and Accurate Ridge Regression via Expectation Maximization
We present a novel method for tuning the regularization hyper-parameter, $λ$, 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 $λ$, 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 $λ$ 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 $λ$ 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 $λ$ (in the regime $q, p \in O(\sqrt{n})$ where $q$ is the number of regression targets).
Paper
Similar papers
Peer review
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.
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. The 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: 1. 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. 2. 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. Best, fBDz William T. Stephenson, Zachary Frangella, Madeleine Udell, Tamara Broderick. Can we globally optimize cross-validation loss? Quasiconvexity in ridge regression. NeurIPS 2021. --------------------------------------------------------------- I 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). My 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. Overall, 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?) I had a few other more contained comments: - I don't think $q$ was defined mathematically in the paper. - 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?) - 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? - 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. - 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! - 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. -In equation (12), it should be clarified that the $\tau^2 (p - n)_+$ is outside the summation. - 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. - 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? I have a few other
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.
Rating
6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
Confidence
4: You are confident in your assessment, but not absolutely certain. It is unlikely, but not impossible, that you did not understand some parts of the submission or that you are unfamiliar with some pieces of related work.
Soundness
4 excellent
Presentation
4 excellent
Contribution
3 good
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. I don't think there is any potential negative societal impact of this work.
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. The 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.
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. Quality: All the numerical experiments seem to be certainly conducted. Numerical codes are also distributed. In this sense, the quality is high. Clarity: 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.
Weaknesses
- The basic idea of the research is not surprising. - 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.).
Questions
- What $k$ and $q$? Please explain about them when they first appear in the text.
Rating
5: Borderline accept: Technically solid paper where reasons to accept outweigh reasons to reject, e.g., limited evaluation. Please use sparingly.
Confidence
3: You are fairly confident in your assessment. It is possible that you did not understand some parts of the submission or that you are unfamiliar with some pieces of related work. Math/other details were not carefully checked.
Soundness
4 excellent
Presentation
3 good
Contribution
2 fair
Limitations
The authors address their research limitations well. I think there is no concern about the potential societal impact.
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.
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. Some 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. It 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. Some of the language and explanation could be tidied up slightly. Some examples: - "that minimizes (3.” - typo? - 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. - “Described in 5”? “Described in Section 5” maybe? - I assume the third order models include linear/quadratic terms and similar, but it’s slightly unclear.
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? What 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? cv.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? Are 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)? If 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? Are 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. LOOCV 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? How 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?
Rating
7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
Confidence
4: You are confident in your assessment, but not absolutely certain. It is unlikely, but not impossible, that you did not understand some parts of the submission or that you are unfamiliar with some pieces of related work.
Soundness
3 good
Presentation
3 good
Contribution
3 good
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.
Summary
This paper studies the comparison of cross validation (CV) and expectation maximization (EM) on tuning the regularization hyper-parameter $\lambda$. The authors explained the running time of these two methods and showed that EM is faster than CV. Then, 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.
Strengths
The proposed EM approach for tuning $\lambda$ is quite interesting. I have not seen this approach before.
Weaknesses
The result is mostly experimental. Most of the theoretical analyses are on the running time. The authors may want to provide some theoretical explanation on why these two methods end up having comparable predictive performances.
Questions
See weaknesses
Rating
5: Borderline accept: Technically solid paper where reasons to accept outweigh reasons to reject, e.g., limited evaluation. Please use sparingly.
Confidence
2: You are willing to defend your assessment, but it is quite likely that you did not understand the central parts of the submission or that you are unfamiliar with some pieces of related work. Math/other details were not carefully checked.
Soundness
3 good
Presentation
3 good
Contribution
2 fair
Limitations
na
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 and 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.
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?
Questions
No questions at this point.
Rating
7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
Confidence
4: You are confident in your assessment, but not absolutely certain. It is unlikely, but not impossible, that you did not understand some parts of the submission or that you are unfamiliar with some pieces of related work.
Soundness
3 good
Presentation
3 good
Contribution
3 good
Limitations
See my remark on missing discussions and point that require further elaboration.
Thank you for the clarification and global response.
I think the results with your clarifications (especially the missingness) now make more sense after your comments. Some of my questions were trying to formalise which elements of the specified method were modelling decisions and which are inferential decisions: the use of hyperpriors in a Bayesian formalism implies hierarchical relationships between parameters that don't have clear correspondence in the regularized frequentist framework of cv.glmnet. It would be instructive for interpretation and future methods development to demonstrate a nuanced understand of which features of the results are a consequence of which features in that respect, e.g. would a Bayesian hierarchical formulation with a Monte Carlo algorithm underpinning the inference provide qualitatively similar results when run to convergence? This may be a minor quibble that should not seriously compromise the paper, but I would appreciate being explored a little more. I think it would further illuminate the discussion surrounding non-i.i.d. data, or correlated or non-Gaussian-distributed betas.
Thanks for the response. I will take this into consideration during the reviewer-AC discussion.
Thank you for the reply! I do see the strength of this kind of method over LOOCV when leaving out a single datapoint is hugely impactful on the solution -- I actually think it could be really interesting to explore that regime further (is there a more general setting than just $X = I_p$ when LOOCV fails and this method still works great? Like when $X \approx I_p$?). And I totally agree with the emphasis of the difficulty in picking the grid of parameter values for LOOCV. What I'm still hung up on is why we shouldn't see gradient descent on the LOOCV objective as a more standard and equally strong approach to selecting $\lambda$. In the above responses, it seems the authors have three points: 1. The LOOCV objective may not be unimodal, whereas the EM objective (the posterior distribution) converges to something unimodal (a normal distribution by the Bernstein-von Mises theorem). 2. If $m = min(n,p)$, then the EM algorithm takes $O(m)$, whereas gradient descent on LOOCV takes $O(nm)$ time. 3. The EM algorithm is proved to be minimax in the paper, whereas it is not known that the optimal LOOCV parameter is minimax. I don't currently think I agree with all of these arguments; I wrote my concerns below with the same numbering: 1. The Berstein-von Mises theorem says that a posterior distribution will converge in total variation distance to a normal distribution as the amount of data goes to infinity (under regularity conditions that are satisfied here). The normal distribution is certainly log-concave (so, the EM objective on it will be convex). But does this mean that the posterior $P_n$ is convex for any $n$? That is, is the following lemma true? **Lemma (maybe)** Suppose $P$ is a log-concave distribution. Further suppose we have a series of distributions $P_n$ indexed by $n$ such that $P_n \to P$ as $n \to \infty$, where the convergence is in total variation distance. Then there exists some $n_0$ such that $P_n$ is log-concave for all $n \geq n_0$ If this lemma isn't true, then I don't see why the convergence of the posterior to a normal distribution helps us with optimization. Do the authors have a reference or proof of this lemma? 2. I agree, the proposed algorithm is faster than LOOCV via gradient descent per iteration. But the number of required iterations of the proposed algorithm is not known to be $o(n)$, so it's hard to compare them theoretically. 3. There is work showing that the risk of the estimator tuned via LOOCV is asymptotically equivalent to that of the optimally tuned ridge estimator. For example, see Theorem 7 of Hastie et al. (2022). Maybe I'm misunderstanding the authors' point about the benefits of minimaxity, but this seems like the best possible guarantee we'd hope for about LOOCV. So, overall, I don't see a theoretical improvement over LOOCV via gradient descent. I also don't see a methodological improvement over LOOCV via gradient descent, since gradient descent does not require a grid, and both EM and gradient descent will require some stopping criteria. There could be an empirical improvement over gradient descent, but I think that has to be demonstrated in the paper. I would find it pretty convincing if the authors could provide a detailed experiment(s) showing that their EM algorithm is significantly faster than LOOCV without any real loss in accuracy. Ideally the "significantly faster" would be at least a couple orders of magnitude to help confirm the conjectures about runtimes on lines 206-211 (i.e. if those conjectures are true, then EM's advantage should become huge for large $n$). Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani. Surprises in High-Dimensional Ridgeless Least Squares Interpolation. Annals of Statistics. 2022.
Thank you for these further thoughtful comments. Regarding point 1, you are correct that the asymptotic normality merely gives us consistency but does not, or at least not obviously imply unimodality for any finite $n$. Motivated by your comment, we have now added a further theorem that establishes that we have unimodality for sufficiently large $n$ when restricting the domain of the posterior to $\tau^2 \geq \epsilon > 0$. Together with the asymptotic concentration of the posterior around the optimal $\boldsymbol \beta_0$ this implies that the EM algorithm, outside of the case $\boldsymbol\beta_0 = {\bf 0}$ identifies the corresponding mode for finite $n$ when initialised with large enough $\tau$. See more details and a proof [in this comment to the global rebuttal](https://openreview.net/forum?id=Ih2yL7o2Gq¬eId=tO1IFaACBl). Regarding point 3, note that the minimaxity criterion is required to hold for all $n$, rather than only asymptotically. Therefore, the asymptotic optimality of the LOOCV approach does not offer any indication regarding its finite sample minimaxity. While our proof does not offer guarantees for arbitrary designs, it does prove minimaxity for the special case of multiple means for all sample sizes $n$ and can be extended to orthogonal designs. Regarding point 2, it is true that we currently do not have a bound on the required number of iterations of the EM algorithm and also not on a gradient descent adoption for LOOCV. However, for the former the empirical evidence strongly indicates a number of iterations approaching a small constant for large $n$, and this presumably would be the ideal behaviour of gradient descent based LOOCV optimisation with well-calibrated step size function, which is yet to be defined. Of course, as discussed, there are problems with LOOCV outside of the employed optimisation method, in particular in sparse covariate settings where LOOCV becomes uninformative for $\lambda$ selection. Note that the multiple means problem (${\bf X}={\bf I}_n$) is only the most extreme such setting, and we have indeed provided in Figure 1 / Section 6.1 a relaxed scenario where ${\bf X}$ is merely sparse but not identity. We concur that it is an important open direction to further characterise this failure of LOOCV. Overall, we fully agree that there are several interesting questions to be explored around the LOOCV objective itself as well as around its optimisation, e.g., a suitable adaption of gradient descent. However, note that our work already uses as a theoretical and practical benchmark an improvement of the state-of-the-art optimisation method for LOOCV. Even the very recent work referenced by you (Hastie et al. 2022), which does include a shortcut method for LOOCV ridge optimisation, appears to be unaware of the further SVD-based speedup that we derive. In conjunction with deriving a novel efficient algorithm for ridge regression parameter estimation, providing solid evidence of its superiority over said improvement of the LOOCV approach in certain regimes, and, now, stimulated by your and other reviewers’ input, also giving a partial theoretical characterisation of the estimator, is that not sufficient for publishing a novel approach to a problem of general interest?
Thank you for the clarification. I have also read other reviewers' comments and found that the proposed method tends to be highly rated, and I admit the proposed method is practically useful. But still, I cannot believe the proposed method is so novel. The Bayesian ridge regression is common, the EM algorithm is common in Bayesian, and SVD utilization at preprocessing is common... Combining these three common pieces makes non-trivial? Maybe yes. I have tried to find some reference writing about this but could not find it. Still I don't know if this is appropriate to be accepted as a theoretical result in NeurIPS, but since I could not provide any direct evidence against the novelty and the proposed method's practicality is clear in terms of the computational cost, I raise my score up to the positive side.
The asymptotic result given in our rebuttal does imply consistency of the posterior, but, as pointed out by reviewer fBDz, it does not, or at least not directly, imply simple posteriors for any finite $n$. Here we further clarify the question regarding the behavior for finite $n$ with an additional result that says that for every $\epsilon > 0$ the posterior, when restricted to $\tau^2 \geq \epsilon$, is unimodal for sufficiently large $n$. This is under the technical assumption that the smallest eigenvalue of ${\bf X'X}$ grows at least linearly with $n$, which is, e.g., satisfied with high probability when the smallest marginal covariate variance is bounded away from $0$. Together with the posterior convergence to normality around the optimal $\boldsymbol \beta_0$, this means that, in the non-degenerate case that $\boldsymbol \beta_0 \neq 0$, there is $\tau^2 = \epsilon > 0$, such that for large enough but finite $n$ the posterior concentrates around $(\boldsymbol \beta_0, \tau^2)$, and thus the corresponding mode is found by EM if initialized with large enough $\tau^2$. Below is the formal statement with proof. **Theorem.** Let $\epsilon > 0$ and $\gamma$ the smallest eigenvalue of ${\bf X'X}/n$. If $\gamma > 0$ and $n > 4/(\gamma\epsilon)$ then the joint posterior $p(\boldsymbol \beta, \sigma^2, \tau^2 | \bf y)$ has exactly one mode with $\tau^2 \geq \epsilon$. *Proof.* We prove the claim by showing that, for sufficiently large $n$, a continuous injective reparameterization of the negative log joint posterior is convex when restricted to $\tau^2 \geq \epsilon$. This is sufficient, since unimodality is preserved by strictly monotone transformations and continuous injective reparameterizations. Specifically, for the presented hierarchical model in (4), the negative log joint posterior up to an additive constant is $$ \frac{n+p+2}{2} \log \sigma^2 + \frac{1}{2\sigma^2}\\|{\bf y} - {\bf X} \boldsymbol \beta\\|^2 + \frac{p+1}{2} \log \tau^2 + \frac{\\|\boldsymbol \beta\\|^2}{2\sigma^2\tau^2} + \log(1 + \tau^2) $$ and reparameterising with $\phi=\beta/\sigma$, $\rho=1/\sigma$ and $\chi = 1/\tau$ and reorganising terms yields $$ -(n+p+2) \log \rho - (p+1) \log \chi + \log(1+\chi^{-2}) + \frac{1}{2}\\|\rho {\bf y}- {\bf X} \boldsymbol \phi\\|^2 + \frac{\\|\chi \boldsymbol \phi\\|^2}{2} . $$ The first three terms can easily be checked to be convex via a second derivative test. For the last two terms, the combined Hessian is of block form $[ {\bf A, B; B', C}]$ with ${\bf A = X'X}+\chi^2 {\bf I}$, ${\bf B}=2\chi \boldsymbol \phi$, and ${\bf C}=\\|\boldsymbol \phi\\|^2$. Symmetric matrices of this form are positive definite if ${\bf A}$ and its Schur complement $$ {\bf C-B'A^{-1}B} = \\|\boldsymbol \phi\\|^2 - 4\boldsymbol \phi'({\bf X'X}/\chi^2+I)^{-1}\boldsymbol \phi $$ are positive definite. Clearly, ${\bf A}$ is positive definite. Moreover, for $n > 4/(\epsilon \gamma)$, we have $$ \begin{aligned} \boldsymbol \phi'({\bf X'X}/\chi^2+I)^{-1}\boldsymbol \phi &= \boldsymbol \phi'({\bf V} \boldsymbol \Sigma' \boldsymbol \Sigma {\bf V}'/\chi^2+I)^{-1}\boldsymbol \phi \\\\ &= \sum_{i=1}^p \phi_i^2 \chi^2 / (s^2_i+\chi^2) \\\\ &< \\|\boldsymbol \phi\\|^2 / (\epsilon n\gamma) \\\\ &< \\|\boldsymbol \phi\\|^2 / 4 \end{aligned} $$ implying the required positivity of the Schur complement.
Decision
Accept (poster)