scikit-learn-contrib / skglm

Fast and modular sklearn replacement for generalized linear models
http://contrib.scikit-learn.org/skglm
BSD 3-Clause "New" or "Revised" License
156 stars 30 forks source link

Unable to fit estimator with Cox datafit and SLOPE penalty #177

Closed tobias-zehnder closed 1 year ago

tobias-zehnder commented 1 year ago

SLOPE is a subclass of BasePenalty, for which is_penalized() is defined as an abstract method. However, SLOPE does not implement is_penalized(). This is not recognized during class instantiation because BasePenalty is not properly implemented as a subclass of ABC and therefore does not enforce the @abstractmethod functionality.

mathurinm commented 1 year ago

Hi @tobias-zehnder, thanks for the catch!

The obvious fix is to implement is_penalized for the SLOPE class, but actually, is_penalized is used only for working set solvers, and there is currently no WS methods suited for SLOPE as the latter is not separable, so this method would not be used in practice

Did you run into this while running some code ? A better fix would be to tackle #101 or something of the like to check that a given solver can support a combination of datafit and penalty.

tobias-zehnder commented 1 year ago

Hi @mathurinm, thanks for the quick reply. Yes, I am trying to run a combination of Cox datafit and SLOPE penalty using the GeneralizedLinearEstimator.

mathurinm commented 1 year ago

Can you share the full snippet? (you can omit the data)

tobias-zehnder commented 1 year ago
# Import modules
import numpy as np
from skglm.estimators import GeneralizedLinearEstimator
from skglm.penalties import SLOPE
from skglm.datafits import Cox

# Define random data
n_samples = 100
n_features = 500
X = np.random.rand(n_samples, n_features)
y = np.array([np.random.rand(n_samples),
              np.random.choice([1,0], n_samples)]).T
alphas = sorted(np.random.rand(n_features), reverse=True)

# Initialize the Cox datafit object
cox_datafit = Cox(use_efron=False)
cox_datafit.initialize(X, y)

# Define the SLOPE penalty
slope_penalty = SLOPE(alphas=alphas)

# Create the GeneralizedLinearEstimator with Cox datafit and SLOPE penalty
estimator = GeneralizedLinearEstimator(datafit=cox_datafit, penalty=slope_penalty)

# Fit the model
model = estimator.fit(X, y)
qklopfenstein-owkin commented 1 year ago

Hello @mathurinm and @tobias-zehnder ! I think the combination Cox datafit + SLOPE is not supported yet. Penalized cox model are solved with Prox Newton but SLOPE does not implement the is_penalized nor subdiff_distance and we cannot use FISTA with Cox datafit because the global lipschitz constant is not computable

Badr-MOUFAD commented 1 year ago

That's right @qklopfenstein-owkin, thanks for flagging that!

To enable combining Cox + SLOPE, I suggest two solutions

qklopfenstein-owkin commented 1 year ago

I can do a quick PR for the upper bound. Do you have a ref for the calculation?

Badr-MOUFAD commented 1 year ago

I shall update the equations in Cox tutorial to account for that after equation 6.

mathurinm commented 1 year ago
  • Implement global_lipschitz for Cox datafit which I believe is upper-bounded by $\frac1n \Vert B^\top s \Vert_\infty \Vert X \Vert_2^2$

Yes this looks correct, since $B_{jj} = 1$ $e^{uj} / B e^u \leq 1$. Also $\Vert B^\top s \Vert\infty$ should be equal to $\displaystyle \sum_{i=1}^n s_i$ since $s$ is positive and $B$ has a least one column full of ones (corresponding to the largest $y_j$).