ENH add Huber datafit

[EnLAI111]

When calling the function construct_grad(X, y, w, Xw, datafit, ws), I get the follow error.

This error may have been caused by the following argument(s):
- argument 4: Cannot determine Numba type of <class '__main__.Huber'>

So maybe there are bugs in the class Huber somewhere, but I can't figure out where.

[PABannier]

Hello @EnLAI111 , you need to decorate your Huber datafit with @jitclass decorator and specify the data type of your class attributes. Look at Quadratic for an example!

[mathurinm]

Thanks a lot for the PR @EnLAI111 ! I did a first pass of comments, don't hesitate to ping me if something is not explicit enough or if you need feedback.

[mathurinm]

I think the lipschitz constant should just be norm(X_j) ** 2?

[QB3]

I would have said norm(X_j) ** 2 / delta.

This makes me think on how do we want to solve this optimization problem?
If we want to use algorithms like in http://proceedings.mlr.press/v84/massias18a/massias18a.pdf and in https://arxiv.org/pdf/1902.02509.pdf, I think we will need a custom _cd_epoch function

[mathurinm]

Sure for delta ? https://github.com/scikit-learn-contrib/skglm/pull/14/files#diff-518236c0559dd6839714e8c437731d42952421ed7fed1319945a7d9bbe9f315eR20

There is no optimization over delta so far so no need for this solver, but if we need to optimize over this variable a dedicated solver will be needed

[QB3]

ok I see, there is no delta for this parametrization!

[EnLAI111]

Thanks a lot for the PR @EnLAI111 ! I did a first pass of comments, don't hesitate to ping me if something is not explicit enough or if you need feedback.

Thanks for the comments and suggestion! I have fixed the problems. Now with WeightedL1 there is no more posted errors of if __name__ == '__main__' execution.

[tomMoral]

Hello @EnLAI111 , you need to decorate your Huber datafit with @jitclass decorator and specify the data type of your class attributes. Look at Quadratic for an example!

Note that this is not explicit in the doc. When reading it, I though it would be automatized by some jit_factaory you have (which was surprising indeed ^^)

Adding the spec_quadratic in the example and explicitely mentionning the need for a jitclass (with a check in the subsequent class) would help I think.

[PABannier]

To finalize this PR, unit tests should be added. Initially I wanted to compare it to sklearn's HuberRegressor but it turns out sklearn optimizes an additional parameter sigma. Any ideas on the unit tests we could write? @mathurinm @QB3

[Klopfe]

To finalize this PR, unit tests should be added. Initially I wanted to compare it to sklearn's HuberRegressor but it turns out sklearn optimizes an additional parameter sigma. Any ideas on the unit tests we could write? @mathurinm @QB3

if we just want to test the convergence, we could use cvxpy to solve the optimization problem and check if the betas are the same. What do you think?

[mathurinm]

@Klopfe using cvxpy would introduce heavy dependency, instead we can fit a sklearn Huber regressor, find its scale, and then fit our model with delta = clf.scale_ * clf.epsilon.

I have pushed a script that does it. Beware that sklearn uses squared L2 regularization, and that the penalty is not scaled by 1/2 (their docstring is not very clear). I had to use more samples than features otherwise it fits perfectly and scale goes to 0

We don't have a way to handle unpenalized problems so I have used WeightedL1 with 0 weights. For future works...

[tomMoral]

another solution could be to use scipy.optimize to solve the problem ?

[mathurinm]

It's what HuberRegressor does already, so why rewrite this code ?

[mathurinm]

Did my pass, merge when green if LGTY @QB3 @PABannier

[mathurinm]

Thanks @EnLAI111 !