implemented lipschitz for dense case, support for AndersonCD and adressed feedback comments

Author: floriankozikowski

examples/plot_smooth_quantile.py

@@ -16,43 +16,18 @@ def pinball_loss(residuals, quantile):
     return np.mean(residuals * (quantile - (residuals < 0)))
 
 
-def create_data(n_samples=1000, n_features=10, noise=0.1):
-    X, y = make_regression(n_samples=n_samples, n_features=n_features, noise=noise)
-    return X, y
-
-
-def plot_quantile_huber():
-    quantiles = [0.1, 0.5, 0.9]
-    delta = 0.5
-    residuals = np.linspace(-3, 3, 500)
-    _, (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]
-        grad = [qh._grad_scalar(r) for r in residuals]
-        ax1.plot(residuals, loss, label=f"τ={tau}")
-        ax2.plot(residuals, grad, label=f"τ={tau}")
-    ax1.set_title("QuantileHuber Loss")
-    ax1.set_xlabel("Residual")
-    ax1.set_ylabel("Loss")
-    ax1.legend()
-    ax2.set_title("QuantileHuber Gradient")
-    ax2.set_xlabel("Residual")
-    ax2.set_ylabel("Gradient")
-    ax2.legend()
-    plt.tight_layout()
-    plt.show()
-
-
 if __name__ == "__main__":
-    X, y = create_data()
+    X, y = make_regression(n_samples=10000, n_features=1000, 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, y)
+    sk.fit(X_c, y_c)
+    sk_pred = sk.predict(X_c) + q_tau
     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)
 
@@ -61,13 +36,14 @@ def plot_quantile_huber():
         quantile=tau,
         alpha=0.1,
         delta_init=0.5,
-        delta_final=0.05,
+        delta_final=0.01,
         n_deltas=5,
+        solver="AndersonCD",
         verbose=True
     )
-    qh.fit(X, y)
+    qh.fit(X_c, y_c)
     qh_time = time.time() - start
-    qh_pred = qh.predict(X)
+    qh_pred = qh.predict(X_c) + q_tau
     qh_cov = np.mean(y <= qh_pred)
     qh_pinball = pinball_loss(y - qh_pred, tau)
 
@@ -79,4 +55,24 @@ def plot_quantile_huber():
     print(f"{'QuantileHuber':<12} {tau:<4} {qh_cov:<8.3f} {qh_time:<6.3f} "
           f"{qh_pinball:<8.4f}")
 
-    plot_quantile_huber()
+
+quantiles = [0.1, 0.5, 0.9]
+delta = 0.5
+residuals = np.linspace(-3, 3, 500)
+_, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
+for tau in quantiles:
+    qh = QuantileHuber(quantile=tau, delta=delta)
+    loss = [qh._loss_sample(r) for r in residuals]
+    grad = [qh._grad_per_sample(r) for r in residuals]
+    ax1.plot(residuals, loss, label=f"τ={tau}")
+    ax2.plot(residuals, grad, label=f"τ={tau}")
+ax1.set_title("QuantileHuber Loss")
+ax1.set_xlabel("Residual")
+ax1.set_ylabel("Loss")
+ax1.legend()
+ax2.set_title("QuantileHuber Gradient")
+ax2.set_xlabel("Residual")
+ax2.set_ylabel("Gradient")
+ax2.legend()
+plt.tight_layout()
+plt.show()

skglm/experimental/quantile_huber.py

@@ -1,14 +1,14 @@
 import numpy as np
+from numpy.linalg import norm
 from numba import float64
-from skglm.datafits.single_task import Huber
+from skglm.datafits.base import BaseDatafit
 from sklearn.base import BaseEstimator, RegressorMixin
-from sklearn.utils.validation import check_X_y, check_array
-from skglm.solvers import FISTA
+from skglm.solvers import FISTA, AndersonCD
 from skglm.penalties import L1
 from skglm.estimators import GeneralizedLinearEstimator
 
 
-class QuantileHuber(Huber):
+class QuantileHuber(BaseDatafit):
     r"""Quantile Huber loss for quantile regression.
 
     Implements the smoothed pinball loss:
@@ -51,11 +51,11 @@ def value(self, y, w, Xw):
         res = 0.0
         for i in range(n_samples):
             residual = y[i] - Xw[i]
-            res += self._loss_scalar(residual)
+            res += self._loss_sample(residual)
         return res / n_samples
 
-    def _loss_scalar(self, residual):
-        """Calculate loss for a single residual."""
+    def _loss_sample(self, residual):
+        """Calculate loss for a single sample."""
         tau = self.quantile
         delta = self.delta
         r = residual
@@ -79,11 +79,11 @@ def gradient_scalar(self, X, y, w, Xw, j):
         grad_j = 0.0
         for i in range(n_samples):
             residual = y[i] - Xw[i]
-            grad_j += -X[i, j] * self._grad_scalar(residual)
+            grad_j += -X[i, j] * self._grad_per_sample(residual)
         return grad_j / n_samples
 
-    def _grad_scalar(self, residual):
-        """Calculate gradient for a single residual."""
+    def _grad_per_sample(self, residual):
+        """Calculate gradient for a single sample."""
         tau = self.quantile
         delta = self.delta
         r = residual
@@ -101,12 +101,26 @@ def _grad_scalar(self, residual):
             # Lower linear tail: r <= -delta
             return tau - 1
 
+    def get_lipschitz(self, X, y):
+        n_features = X.shape[1]
+
+        lipschitz = np.zeros(n_features, dtype=X.dtype)
+        c = max(self.quantile, 1 - self.quantile) / self.delta
+        for j in range(n_features):
+            lipschitz[j] = c * (X[:, j] ** 2).sum() / len(y)
+
+        return lipschitz
+
+    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)
+
 
 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):
+                 n_deltas=10, max_iter=1000, tol=1e-4, verbose=False, solver="FISTA"):
         self.quantile = quantile
         self.alpha = alpha
         self.delta_init = delta_init
@@ -115,10 +129,10 @@ def __init__(self, quantile=0.5, alpha=0.1, delta_init=1.0, delta_final=1e-3,
         self.max_iter = max_iter
         self.tol = tol
         self.verbose = verbose
+        self.solver = solver
 
     def fit(self, X, y):
         """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)
 
@@ -127,19 +141,26 @@ def fit(self, X, y):
                 f"Progressive smoothing: delta {self.delta_init:.3f} --> "
                 f"{self.delta_final:.3f} in {self.n_deltas} steps")
 
-        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)
+        datafit = QuantileHuber(quantile=self.quantile, delta=self.delta_init)
+        penalty = L1(alpha=self.alpha)
+        # Solver selection
+        if isinstance(self.solver, str):
+            if self.solver == "FISTA":
+                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)
+            else:
+                raise ValueError(f"Unknown solver: {self.solver}")
+        else:
+            solver = self.solver
 
-            est = GeneralizedLinearEstimator(
-                datafit=datafit,
-                penalty=penalty,
-                solver=solver
-            )
+        est = GeneralizedLinearEstimator(
+            datafit=datafit, penalty=penalty, solver=solver)
 
-            if i > 0:
-                est.coef_ = w.copy()
+        for i, delta in enumerate(deltas):
+            datafit.delta = float(delta)
 
             est.fit(X, y)
             w = est.coef_.copy()
@@ -151,13 +172,16 @@ def fit(self, X, y):
 
                 print(
                     f"  Stage {i+1:2d}: delta={delta:.4f}, "
-                    f"coverage={coverage:.3f}, pinball_loss={pinball_loss:.6f}")
+                    f"coverage={coverage:.3f}, pinball_loss={pinball_loss:.6f}, "
+                    f"n_iter={est.n_iter_}"
+                )
 
-        self.coef_ = w
+        self.est = est
 
         return self
 
     def predict(self, X):
         """Predict using the fitted model."""
-        X = check_array(X)
-        return X @ self.coef_
+        if not hasattr(self, "est"):
+            raise ValueError("Call 'fit' before 'predict'.")
+        return self.est.predict(X)