As reported in #13923 the currently used scipy.linalg.lstsq can be significantly slower than Ridge(solver="cholesky", alpha=0) for tall dense X.

We should also call scipy.linalg.lstsq with check_finite=False as we already do input validation in fit.

Also the scipy.linalg.lstsq function has an optional lapack_driver that accepts the following options: 'gelsd', 'gelsy', 'gelss'. The default is'gelsd'. Maybe we should expose the others in our API and benchmark them to see if it would make sense to change the default.

Related issue for Ridge: #14269