add intercept method (only works for AndersonCD so far)

Author: floriankozikowski

examples/plot_smooth_quantile.py

@@ -8,24 +8,19 @@
 import matplotlib.pyplot as plt
 from sklearn.datasets import make_regression
 
-# TODO: no intercept handling yet
-
 
 def pinball_loss(residuals, quantile):
     """True pinball loss."""
     return np.mean(residuals * (quantile - (residuals < 0)))
 
 
-X, y = make_regression(n_samples=10000, n_features=1000, noise=0.1, random_state=0)
+X, y = make_regression(n_samples=1000, n_features=10, noise=0.1, random_state=0)
 tau = 0.8
-X_c = X - X.mean(axis=0)
-q_tau = np.quantile(y, tau)
-y_c = y - q_tau
 
 start = time.time()
-sk = QuantileRegressor(quantile=tau, alpha=0.1, fit_intercept=False)
-sk.fit(X_c, y_c)
-sk_pred = sk.predict(X_c) + q_tau
+sk = QuantileRegressor(quantile=tau, alpha=0.1, fit_intercept=True)
+sk.fit(X, y)
+sk_pred = sk.predict(X)
 sk_time = time.time() - start
 sk_cov = np.mean(y <= sk_pred)
 sk_pinball = pinball_loss(y - sk_pred, tau)
@@ -38,11 +33,12 @@ def pinball_loss(residuals, quantile):
     delta_final=0.01,
     n_deltas=5,
     solver="AndersonCD",
-    verbose=True
+    verbose=True,
+    fit_intercept=True
 )
-qh.fit(X_c, y_c)
+qh.fit(X, y)
 qh_time = time.time() - start
-qh_pred = qh.predict(X_c) + q_tau
+qh_pred = qh.predict(X)
 qh_cov = np.mean(y <= qh_pred)
 qh_pinball = pinball_loss(y - qh_pred, tau)
 

skglm/experimental/quantile_huber.py

@@ -115,12 +115,21 @@ def get_global_lipschitz(self, X, y):
         c = max(self.quantile, 1 - self.quantile) / self.delta
         return c * norm(X, ord=2) ** 2 / len(y)
 
+    def intercept_update_step(self, y, Xw):
+        n_samples = len(y)
+        update = 0.0
+        for i in range(n_samples):
+            residual = y[i] - Xw[i]
+            update -= self._grad_per_sample(residual)
+        return update / n_samples
+
 
 class SmoothQuantileRegressor(BaseEstimator, RegressorMixin):
     """Quantile regression with progressive smoothing."""
 
     def __init__(self, quantile=0.5, alpha=0.1, delta_init=1.0, delta_final=1e-3,
-                 n_deltas=10, max_iter=1000, tol=1e-4, verbose=False, solver="FISTA"):
+                 n_deltas=10, max_iter=1000, tol=1e-4, verbose=False,
+                 solver="AndersonCD", fit_intercept=True):
         self.quantile = quantile
         self.alpha = alpha
         self.delta_init = delta_init
@@ -130,6 +139,7 @@ def __init__(self, quantile=0.5, alpha=0.1, delta_init=1.0, delta_final=1e-3,
         self.tol = tol
         self.verbose = verbose
         self.solver = solver
+        self.fit_intercept = fit_intercept
 
     def fit(self, X, y):
         """Fit using progressive smoothing: delta_init --> delta_final."""
@@ -146,11 +156,18 @@ def fit(self, X, y):
         # Solver selection
         if isinstance(self.solver, str):
             if self.solver == "FISTA":
+                if self.fit_intercept:
+                    import warnings
+                    warnings.warn(
+                        "FISTA solver does not support intercept. "
+                        "Falling back to fit_intercept=False."
+                    )
+                    self.fit_intercept = False
                 solver = FISTA(max_iter=self.max_iter, tol=self.tol)
                 solver.warm_start = True
             elif self.solver == "AndersonCD":
                 solver = AndersonCD(max_iter=self.max_iter, tol=self.tol,
-                                    warm_start=True, fit_intercept=False)
+                                    warm_start=True, fit_intercept=self.fit_intercept)
             else:
                 raise ValueError(f"Unknown solver: {self.solver}")
         else:
@@ -167,6 +184,8 @@ def fit(self, X, y):
 
             if self.verbose:
                 residuals = y - X @ w
+                if self.fit_intercept:
+                    residuals -= est.intercept_
                 coverage = np.mean(residuals <= 0)
                 pinball_loss = np.mean(residuals * (self.quantile - (residuals < 0)))