Newton regularize

our `optimize._fit_newton` is very plain and fails if the Hessian is not invertible

see #953 for the suggestion to use pinv

This PR instead introduces a Ridge regularization, i.e. add `lambda * eye`

I'm not sure how lambda should be set.
The current lambda is at least 1e-10, which I needed to avoid the failing cases in #1845. Making it only depend on median or max of the diagonal elements didn't work in that case.

Problem in emplike #1845 is that _fit_newton is hardcoded as the only optimizer used in the inner optimization.