add progressive smoothing

Author: floriankozikowski

examples/plot_smooth_quantile.py

@@ -4,11 +4,11 @@
 import numpy as np
 import time
 from sklearn.linear_model import QuantileRegressor
-from skglm.experimental.quantile_huber import QuantileHuber, SimpleQuantileRegressor
+from skglm.experimental.quantile_huber import QuantileHuber, SmoothQuantileRegressor
 import matplotlib.pyplot as plt
 from sklearn.datasets import make_regression
 
-# TODO: no smoothing and no intercept handling yet
+# TODO: no intercept handling yet
 
 
 def pinball_loss(residuals, quantile):
@@ -25,7 +25,7 @@ def plot_quantile_huber():
     quantiles = [0.1, 0.5, 0.9]
     delta = 0.5
     residuals = np.linspace(-3, 3, 500)
-    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
+    _, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
     for tau in quantiles:
         qh = QuantileHuber(quantile=tau, delta=delta)
         loss = [qh._loss_scalar(r) for r in residuals]
@@ -49,15 +49,22 @@ def plot_quantile_huber():
     tau = 0.8
 
     start = time.time()
-    sk = QuantileRegressor(quantile=tau, alpha=0.001, fit_intercept=False)
+    sk = QuantileRegressor(quantile=tau, alpha=0.1, fit_intercept=False)
     sk.fit(X, y)
     sk_time = time.time() - start
     sk_pred = sk.predict(X)
     sk_cov = np.mean(y <= sk_pred)
     sk_pinball = pinball_loss(y - sk_pred, tau)
 
     start = time.time()
-    qh = SimpleQuantileRegressor(quantile=tau, alpha=0.001, delta=0.05)
+    qh = SmoothQuantileRegressor(
+        quantile=tau,
+        alpha=0.1,
+        delta_init=0.5,
+        delta_final=0.05,
+        n_deltas=5,
+        verbose=True
+    )
     qh.fit(X, y)
     qh_time = time.time() - start
     qh_pred = qh.predict(X)

skglm/experimental/__init__.py

@@ -2,7 +2,7 @@
 from .sqrt_lasso import SqrtLasso, SqrtQuadratic
 from .pdcd_ws import PDCD_WS
 from .quantile_regression import Pinball
-from .quantile_huber import QuantileHuber, SimpleQuantileRegressor
+from .quantile_huber import QuantileHuber, SmoothQuantileRegressor
 
 __all__ = [
     IterativeReweightedL1,
@@ -11,5 +11,5 @@
     SqrtQuadratic,
     SqrtLasso,
     QuantileHuber,
-    SimpleQuantileRegressor,
+    SmoothQuantileRegressor,
 ]

skglm/experimental/quantile_huber.py

@@ -1,3 +1,4 @@
+import numpy as np
 from numba import float64
 from skglm.datafits.single_task import Huber
 from sklearn.base import BaseEstimator, RegressorMixin
@@ -101,32 +102,58 @@ def _grad_scalar(self, residual):
             return tau - 1
 
 
-class SimpleQuantileRegressor(BaseEstimator, RegressorMixin):
-    """Simple quantile regression without progressive smoothing."""
+class SmoothQuantileRegressor(BaseEstimator, RegressorMixin):
+    """Quantile regression with progressive smoothing."""
 
-    def __init__(self, quantile=0.5, alpha=0.1, delta=0.1, max_iter=1000, tol=1e-4):
+    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):
         self.quantile = quantile
         self.alpha = alpha
-        self.delta = delta
+        self.delta_init = delta_init
+        self.delta_final = delta_final
+        self.n_deltas = n_deltas
         self.max_iter = max_iter
         self.tol = tol
+        self.verbose = verbose
 
     def fit(self, X, y):
-        """Fit using FISTA with fixed delta."""
+        """Fit using progressive smoothing: delta_init --> delta_final."""
         X, y = check_X_y(X, y)
+        w = np.zeros(X.shape[1])
+        deltas = np.geomspace(self.delta_init, self.delta_final, self.n_deltas)
 
-        datafit = QuantileHuber(quantile=self.quantile, delta=self.delta)
-        penalty = L1(alpha=self.alpha)
-        solver = FISTA(max_iter=self.max_iter, tol=self.tol)
+        if self.verbose:
+            print(
+                f"Progressive smoothing: delta {self.delta_init:.3f} --> "
+                f"{self.delta_final:.3f} in {self.n_deltas} steps")
 
-        est = GeneralizedLinearEstimator(
-            datafit=datafit,
-            penalty=penalty,
-            solver=solver
-        )
+        for i, delta in enumerate(deltas):
+            datafit = QuantileHuber(quantile=self.quantile, delta=delta)
+            penalty = L1(alpha=self.alpha)
+            solver = FISTA(max_iter=self.max_iter, tol=self.tol)
 
-        est.fit(X, y)
-        self.coef_ = est.coef_
+            est = GeneralizedLinearEstimator(
+                datafit=datafit,
+                penalty=penalty,
+                solver=solver
+            )
+
+            if i > 0:
+                est.coef_ = w.copy()
+
+            est.fit(X, y)
+            w = est.coef_.copy()
+
+            if self.verbose:
+                residuals = y - X @ w
+                coverage = np.mean(residuals <= 0)
+                pinball_loss = np.mean(residuals * (self.quantile - (residuals < 0)))
+
+                print(
+                    f"  Stage {i+1:2d}: delta={delta:.4f}, "
+                    f"coverage={coverage:.3f}, pinball_loss={pinball_loss:.6f}")
+
+        self.coef_ = w
 
         return self