@agramfort @lorentzenchr @jjerphan I think the clarification of the motivations for the centering in _preprocess_data and post-hoc intercept computation in _set_intercept is ready to review. Here is the rendered docs:

• https://output.circle-artifacts.com/output/job/6b1ee83a-3571-4f3b-930f-e819b7edc5e5/artifacts/0/doc/modules/linear_model.html#note-on-the-handling-of-the-intercept-in-linear-regression-models

I still need to think on how to update the ridge tests to be consistent with this for all cases. I probably won't have time to do it today though.

In my lately short style - but implicitly very appreciating —just a few points:

• I think it does not solve BUG unpenalized Ridge does not give minimum norm solution #22947 because the math for centering with n_features>n_samples does not work out.
• These notes only concern the implementation but should not concern a user. Should we really place a solution to a standard exercise of a stats/ml lecture in out user guide?

I think it does not solve #22947 because the math for centering with n_features>n_samples does not work out.

Can you be more precise?

EDIT: indeed there is a problem when I wrote "we recover:". There could be other values of w_0 that cause the sample wise sum equality to hold.

These notes only concern the implementation but should not concern a user. Should we really place a solution to a standard exercise of a stats/ml lecture in out user guide?

I agree. We can move them to another place (once fixed). But I wanted to benefit from the LaTeX rendering. An ascii version would be hard to read in my opinion.

Maybe we could move that to a maintainer-oriented section of the doc.

I have a proof sketch for the under-determined case using the method of Lagrange multipliers (for the minimization problem that does not include w_0 in the computation of the norm). I will try to finalize it on Monday.

Ok I pushed my proof for the underdetermined case here:

https://github.com/ogrisel/minimum-norm-ols/blob/main/minimum-norm-ols-intercept.pdf

/cc @lorentzenchr @agramfort @jjerphan

I will update the restructured text versions of the notes based on your feedback.

Thank you for writing this down, @ogrisel.

To me, the derivation is correct, and we just need an argument showing that $X_c X_c^\top$ is invertible.

(Edit: also since is convex, the minimizer necessarily is the critical point and we directly can proceed with equivalences.)

These notes only concern the implementation but should not concern a user.

The only user impacting thing it that we would now guarantee that LinearRegression(with_intercept=True).fit(X, y).coef_ is the limit of Ridge(with_intercept=True, alpha=alpha).fit(X, y).coef_ when alpha goes to zero which was never made explicit before. But can just write it in the user guide and/or the docstring of Ridge's alpha and LinearRegression's "see also" section.

Should we really place a solution to a standard exercise of a stats/ml lecture in out user guide?

Well, standard lectures do not cover the case where we do not penalize the intercept or include it in the computation of the minimum norm. Making it explicit as I did my linked PDF is boring but not that easy to check, and we got the initial version of the fixture wrong because of this.

To me, the derivation is correct, and we just need an argument showing that
is invertible.

Indeed. The solution would still be valid with the pseudo-inverse instead of the inverse (this is what I have to do in the fixture to make it work with any value of global_random_seed) but I am too lazy to prove it :)

The minimizer of the objective function is a critical point, but the converse is not true generally. Here we proceed with the converse to derive the solution, but we still need an argument to show that there's an equivalence here.

The minimizer of the objective function is a critical point, but the converse is not true generally. Here we proceed with the converse to derive the solution, but we still need an argument to show that there's an equivalence here.

I think the equivalence comes from the Lagrange multiplier theorem:

https://en.wikipedia.org/wiki/Lagrange_multiplier#Statement

Just for the record, I still plan to follow-up on this. Thanks @agramfort for the review. Do you have any opinion on how to document the under-determined case?

Do you think I should port the content of https://github.com/ogrisel/minimum-norm-ols/blob/main/minimum-norm-ols-intercept.pdf into our sphinx doc (assuming you agree with the content)?

@ogrisel I also plan to review this PR. From the glance I took, my main concern is the proposal to not include the intercept in the norm.

@lorentzenchr @agramfort I finally took the time to conduct a numerical study of the impact of the problem formulation:

https://gist.github.com/ogrisel/18bbf02128a3890b03534d52b7fc8bd0#file-minimum_norm_ridge_limit-ipynb

The current formulation (ridge on centered data followed by analytical intercept computation) seems to beneficial for several reasons detailed in the notebooks. The main problems with what I named "type b" are:

• the intercept estimate is numerically discontinuous when changing alpha to very low or large values, at least when keeping the rcond singular value threshold of the underlying SVD-based solver constant.
• the high regularization value of the intercept is wrong, it does not converge to y_train.mean() as expected.
• I don't see how it would be possible to get a continuous intercept when alpha goes to zero if the intercept is included in the computation of the minimum norm solution but not penalized in the ridge problem, even if we had not fixed the rcond parameter.

We could improve the notebook to check that those conclusions hold when changing the solver.

EDIT: I updated the linked notebook to add a ridge solvers based on lstsq instead of pinv and the conclusions still hold (although the lstsq solver seems slightly less stable for a reason I do not understand).

@ogrisel I‘ll have a deeper look at your great analysis. But I need some time. You could even consider to write a paper about it or a blog post.

Thanks. I won't have the time to dig deeper next week, but ideally we could try to understand better why/how the different solvers break on extreme values of the regularization parameter. We could also probably find out analytically why type b ridge (without intercept penalization) cannot reach the type b OLS (with intercept in the norm). We could do that either generally (for any X, y training set) or for some degenerate cases such as a single data point dataset: X = [[1.]], y = [1.0] for instance.

On the shorter term I would like to first finalize the updates of this PR to make sure that all solvers in scikit-learn behave correctly using the centered ridge formulation.

@ogrisel https://gist.github.com/ogrisel/18bbf02128a3890b03534d52b7fc8bd0#file-minimum_norm_ridge_limit-ipynb has a bug. For lstsq, one needs to take the square root of alpha. It then never fails and is faster than pinv, at least for smaller alpha.

I additionally added np.linalg.solve. It fails for alpha=0, but except for this case, it is as stable (well, X here is well behaved) as pinv and faster.

Fix + additional solver in https://github.com/lorentzenchr/notebooks/blob/master/minimum_norm_ridge_limit.ipynb.

My summary
I find the higher test R2 for the centered variant ("a") in addition to the smoother behaviour of the intercept as a function of alpha somehow striking.

✔️ Linting Passed

All linting checks passed. Your pull request is in excellent shape! ☀️

_{Generated for commit: 036552b. Link to the linter CI: here}

@lorentzenchr I updated my copy of the notebook to include your fix + the numpy.linalg.solve solver and further added type "a" and type "b" solutions to the minimum norm OLS and ridge problems with the lsqr solver to the analysis.

For the latter, the results are very similar to what we observe with lstsq: for the type "b" formulation, with observe discontinuities (unstable) solution paths both for low and high values of alpha that are in turn detrimental to the test R2 performance. Type "a" is always stable, whatever the choice of the solver, except numpy.linalg.solve for small alpha values, but this is expected (as per its docstring).

I plan to resume the work on updating the tests in this PR and making them pass for all solvers and all data shapes using the type "a" formulation as reference solution for the minimum norm OLS problem.

Some related parallel fixes for LinearRegression:

• Fix LinearRegression's numerical stability on rank deficient data by setting the cond parameter in the call to scipy.linalg.lstsq #30040
• LinearRegression on sparse matrices is not sample weight consistent #30131