ENH Add FISTA solver

[PABannier]

Closes #89

A few points to discuss:

1. Currently the FISTA solver uses Gram updates (as per Gram-based CD/BCD/FISTA solvers for (group)Lasso when n_samples >> n_features #4 ). Question: do we want to keep it this way? Implement without Gram update? Or have two options?

2. For non-coordinate-wise updates, we run into the issue of not having a prox_vec method in the BasePenalty class. If we want to support a larger class of penalties for FISTA (e.g.: L1, WeightedL1, SLOPE, ...), we need a prox_vec method.

[mathurinm]

• Remove Gram to handle the generic case (+ gram is only suited to quadratics)
• To keep the API simple, can you do a for loop over coordinates to compute the prox, calling prox_1D ? We'll lose a bit of time but the gradient computation should be the dominating cost

[Badr-MOUFAD]

Thanks for the hard work @PABannier!

Below, some minor remarks.

Besides, I have one concern: I don't think it's a good idea to add support for FISTA to all the datafits. We can limit ourselves to one of them just for testing purposes.
Indeed, AndersonCD and ProxNewton are much faster for separable problems. We better keep FISTA for particular cases (e.g SLOPE #92)

WDYT?

[PABannier]

@Badr-MOUFAD totally agree. FISTA would be for a subset of penalties where PN or AndersonCD are not available.

[mathurinm]

that will yield a very crude bound, potentially with a loss or the order of n_features.

Use a few iterations of the power method instead to approximate the lipschitz constant of the sparse matrix (there's also the Lanczos iteration but it's more complicated, let's implement the easy one first)

[Badr-MOUFAD]

I am unable to find the root cause of the problem at the uniitest.

Here is a small script to reproduce

import numpy as np
from  scipy.sparse import random
from skglm.estimators import LinearSVC

n_samples, n_features = 20, 30
X_sparse = random(n_samples, n_features, density=0.5, format='csc', random_state=0)
y = np.ones(n_samples)

LinearSVC(C=1., tol=1e-9).fit(X_sparse, y)

Output (it depends)

Segmentation fault (core dumped)

or

corrupted size vs. prev_size
Aborted (core dumped)

or

python3: malloc.c:3852: _int_malloc: Assertion `chunk_main_arena (fwd)' failed.
Aborted (core dumped)

@mathurinm, @PABannier, any thoughts?

[PABannier]

A segfault is usually thrown by Numba when it can't access something it should (e.g. missing initialization of datafit). Have you tried setting breakpoints at various places of the code to see which line is causing the issue?

[Badr-MOUFAD]

Yes absolutely, I tried that. it breaks down in the initialization of datafit. Yet, I can't figure out why.
What is surprising is that it works for some X sizes,

[PABannier]

Weird, that works for me on this branch. Can you try reinstalling numba and skglm? I've run into similar issues with celer where I had segfaults on my machine that disappeared when I reinstalled the package.

[Badr-MOUFAD]

I found the bug.
In the case of SVC, the design matrix is yXT which has particularly n_rows=n_features instead of n_samples.
This caused an index out of range in spectral_norm.

Thanks @PABannier for your help!

[mathurinm]

Wow, good catch @Badr-MOUFAD

For a more robust design replace n_samples by n_rows_X and pass X.shape[0] explicitly

[mathurinm]

Thanks @PABannier and @Badr-MOUFAD