FISTA with global lipschitz

Author: PABannier

skglm/solvers/fista.py

@@ -33,9 +33,8 @@ def solve(self, X, y, datafit, penalty, w_init=None, Xw_init=None):
         z = w_init.copy() if w_init is not None else np.zeros(n_features)
         Xw = Xw_init.copy() if Xw_init is not None else np.zeros(n_samples)
 
-        # line search?
-        # lipschitz = np.max(datafit.lipschitz)
-        lipschitz = np.linalg.norm(X, ord=2) ** 2 / n_samples
+        # TODO: OR line search
+        lipschitz = datafit.global_lipschitz
 
         for n_iter in range(self.max_iter):
             t_old = t_new

toy_fista.py

@@ -1,9 +1,9 @@
 import numpy as np
 from numpy.linalg import norm
-from skglm.datafits.single_task import Quadratic
+from skglm.datafits.single_task import Quadratic, Logistic
 from skglm.solvers import FISTA
 from skglm.penalties import L1
-from skglm.estimators import Lasso
+from skglm.estimators import SparseLogisticRegression, Lasso
 from skglm.utils import make_correlated_data, compiled_clone
 
 
@@ -19,6 +19,10 @@
 tol = 1e-10
 
 solver = FISTA(max_iter=max_iter, tol=tol, opt_freq=obj_freq, verbose=1)
+
+##############
+# Quadratic
+##############
 datafit = compiled_clone(Quadratic())
 datafit.initialize(X, y)
 penalty = compiled_clone(L1(alpha))
@@ -28,3 +32,20 @@
 clf.fit(X, y)
 
 np.testing.assert_allclose(w, clf.coef_, rtol=1e-5)
+
+##############
+# Logistic
+##############
+y = np.sign(y)
+alpha_max = norm(X.T @ y, ord=np.inf) / (4 * n_samples)
+alpha = alpha_max / 10
+
+datafit = compiled_clone(Logistic())
+datafit.initialize(X, y)
+penalty = compiled_clone(L1(alpha))
+w = solver.solve(X, y, datafit, penalty)
+
+clf = SparseLogisticRegression(alpha=alpha, tol=tol, fit_intercept=False)
+clf.fit(X, y)
+
+np.testing.assert_allclose(w, np.squeeze(clf.coef_), rtol=1e-3)