make basic version without intercept handling and progressive smoothing

Author: floriankozikowski

examples/plot_smooth_quantile.py

@@ -1,75 +1,75 @@
 """
-===========================================
-Smooth Quantile Regression Example
-===========================================
-
+QuantileHuber vs Sklearn
 """
-
 import numpy as np
-import matplotlib.pyplot as plt
 import time
-from sklearn.datasets import make_regression
-from sklearn.preprocessing import StandardScaler
 from sklearn.linear_model import QuantileRegressor
-from skglm.experimental.smooth_quantile_regressor import SmoothQuantileRegressor
-from skglm.experimental.quantile_huber import QuantileHuber
+from skglm.experimental.quantile_huber import QuantileHuber, SimpleQuantileRegressor
+import matplotlib.pyplot as plt
+from sklearn.datasets import make_regression
+
+# TODO: no smoothing and 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=1000, n_features=10, noise=0.1, random_state=42)
-X = StandardScaler().fit_transform(X)
-tau = 0.75
 
-t0 = time.time()
-reg_skglm = SmoothQuantileRegressor(quantile=tau).fit(X, y)
-t1 = time.time()
-reg_sklearn = QuantileRegressor(quantile=tau, alpha=0.1, solver='highs').fit(X, y)
-t2 = time.time()
+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
 
-y_pred_skglm, y_pred_sklearn = reg_skglm.predict(X), reg_sklearn.predict(X)
-coverage_skglm = np.mean(y <= y_pred_skglm)
-coverage_sklearn = np.mean(y <= y_pred_sklearn)
 
-print(f"\nTiming: skglm={t1-t0:.3f}s, sklearn={t2-t1:.3f}s, "
-      f"speedup={(t2-t1)/(t1-t0):.1f}x")
-print(f"Coverage (target {tau}): skglm={coverage_skglm:.3f}, "
-      f"sklearn={coverage_sklearn:.3f}")
-print(f"Non-zero coefs: skglm={np.sum(reg_skglm.coef_ != 0)}, "
-      f"sklearn={np.sum(reg_sklearn.coef_ != 0)}")
+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))
+    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()
 
 
-# Visualizations
-def pinball(y_true, y_pred):
-    diff = y_true - y_pred
-    return np.mean(np.where(diff >= 0, tau * diff, (1 - tau) * -diff))
+if __name__ == "__main__":
+    X, y = create_data()
+    tau = 0.8
 
+    start = time.time()
+    sk = QuantileRegressor(quantile=tau, alpha=0.001, 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)
 
-print(f"Pinball loss: skglm={pinball(y, y_pred_skglm):.4f}, "
-      f"sklearn={pinball(y, y_pred_sklearn):.4f}")
+    start = time.time()
+    qh = SimpleQuantileRegressor(quantile=tau, alpha=0.001, delta=0.05)
+    qh.fit(X, y)
+    qh_time = time.time() - start
+    qh_pred = qh.predict(X)
+    qh_cov = np.mean(y <= qh_pred)
+    qh_pinball = pinball_loss(y - qh_pred, tau)
 
-plt.figure(figsize=(12, 5))
-plt.subplot(121)
-residuals = np.linspace(-2, 2, 1000)
-for delta in [1.0, 0.5, 0.1]:
-    loss = QuantileHuber(quantile=tau, delta=delta)
-    losses = [loss.value(np.array([r]), np.array([[1]]), np.array([0]))
-              for r in residuals]
-    plt.plot(residuals, losses, label=f'δ={delta}')
-plt.plot(residuals, [tau * max(r, 0) + (1 - tau) * max(-r, 0)
-                     for r in residuals], 'k--', label='Pinball')
-plt.axvline(x=0, color='k', linestyle='--', alpha=0.3)
-plt.xlabel('Residual (y - y_pred)')
-plt.ylabel('Loss')
-plt.title('Quantile Huber Loss (τ=0.75)')
-plt.legend()
-plt.grid(True, alpha=0.3)
+    print(f"{'Method':<12} {'Q':<4} {'Coverage':<8} {'Time':<6} "
+          f"{'Pinball':<8}")
+    print("-" * 55)
+    print(f"{'Sklearn':<12} {tau:<4} {sk_cov:<8.3f} {sk_time:<6.3f} "
+          f"{sk_pinball:<8.4f}")
+    print(f"{'QuantileHuber':<12} {tau:<4} {qh_cov:<8.3f} {qh_time:<6.3f} "
+          f"{qh_pinball:<8.4f}")
 
-plt.subplot(122)
-plt.hist(y - y_pred_skglm, bins=50, alpha=0.5, label='skglm')
-plt.hist(y - y_pred_sklearn, bins=50, alpha=0.5, label='sklearn')
-plt.axvline(0, color='k', linestyle='--')
-plt.xlabel('Residual (y - y_pred)')
-plt.ylabel('Count')
-plt.title('Residuals Histogram')
-plt.legend()
-plt.grid(True, alpha=0.3)
-plt.tight_layout()
-plt.show()
+    plot_quantile_huber()

skglm/experimental/__init__.py

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

skglm/experimental/quantile_huber.py

@@ -1,13 +1,16 @@
-import numpy as np
 from numba import float64
 from skglm.datafits.single_task import Huber
-from skglm.utils.sparse_ops import spectral_norm
+from sklearn.base import BaseEstimator, RegressorMixin
+from sklearn.utils.validation import check_X_y, check_array
+from skglm.solvers import FISTA
+from skglm.penalties import L1
+from skglm.estimators import GeneralizedLinearEstimator
 
 
 class QuantileHuber(Huber):
     r"""Quantile Huber loss for quantile regression.
 
-    Implements the smoothed pinball loss with quadratic region:
+    Implements the smoothed pinball loss:
 
     .. math::
 
@@ -25,17 +28,13 @@ class QuantileHuber(Huber):
         Desired quantile level between 0 and 1.
     delta : float, default=1.0
         Width of quadratic region.
-
-    References
-    ----------
-    Chen, C. (2007). A Finite Smoothing Algorithm for Quantile Regression.
-    Journal of Computational and Graphical Statistics, 16(1), 136–164.
-    http://www.jstor.org/stable/27594233
     """
 
     def __init__(self, quantile=0.5, delta=1.0):
         if not 0 < quantile < 1:
             raise ValueError("quantile must be between 0 and 1")
+        if delta <= 0:
+            raise ValueError("delta must be positive")
         self.delta = float(delta)
         self.quantile = float(quantile)
 
@@ -45,80 +44,93 @@ def get_spec(self):
     def params_to_dict(self):
         return dict(delta=self.delta, quantile=self.quantile)
 
-    def _loss_and_grad_scalar(self, residual):
-        """Calculate loss and gradient for a single residual."""
+    def value(self, y, w, Xw):
+        """Compute the quantile Huber loss value."""
+        n_samples = len(y)
+        res = 0.0
+        for i in range(n_samples):
+            residual = y[i] - Xw[i]
+            res += self._loss_scalar(residual)
+        return res / n_samples
+
+    def _loss_scalar(self, residual):
+        """Calculate loss for a single residual."""
         tau = self.quantile
         delta = self.delta
-        abs_r = abs(residual)
-
-        # Quadratic core: |r| ≤ delta
-        if abs_r <= delta:
-            if residual >= 0:
-                # 0 ≤ r ≤ delta
-                loss = tau * residual**2 / (2 * delta)
-                grad = tau * residual / delta
-            else:
-                # -delta ≤ r < 0
-                loss = (1 - tau) * residual**2 / (2 * delta)
-                grad = (1 - tau) * residual / delta
-            return loss, grad
-
-        # Linear tails: |r| > delta
-        if residual > delta:
-            loss = tau * (residual - delta / 2)
-            grad = tau
-            return loss, grad
+        r = residual
+
+        if r >= delta:
+            # Upper linear tail: r >= delta
+            return tau * (r - delta/2)
+        elif r >= 0:
+            # Upper quadratic: 0 <= r < delta
+            return tau * r**2 / (2 * delta)
+        elif r > -delta:
+            # Lower quadratic: -delta < r < 0
+            return (1 - tau) * r**2 / (2 * delta)
         else:
-            loss = (1 - tau) * (-residual - delta / 2)
-            grad = tau - 1
-            return loss, grad
+            # Lower linear tail: r <= -delta
+            return (1 - tau) * (-r - delta/2)
 
-    def value(self, y, w, Xw):
-        """Compute the quantile Huber loss value."""
-        residuals = y - Xw
-        loss = np.zeros_like(residuals)
-        for i, r in enumerate(residuals):
-            loss[i], _ = self._loss_and_grad_scalar(r)
-        return np.mean(loss)
-
-    def raw_grad(self, y, Xw):
-        """Compute gradient of datafit w.r.t Xw."""
-        residuals = y - Xw
-        grad = np.zeros_like(residuals)
-        for i, r in enumerate(residuals):
-            _, grad[i] = self._loss_and_grad_scalar(r)
-        return -grad
-
-    def get_lipschitz(self, X, y):
-        """Compute coordinate-wise Lipschitz constants."""
-        weight = max(self.quantile, 1 - self.quantile)
-        return weight * (X ** 2).sum(axis=0) / (len(y) * self.delta)
-
-    def get_global_lipschitz(self, X, y):
-        """Compute global Lipschitz constant."""
-        weight = max(self.quantile, 1 - self.quantile)
-        return weight * np.linalg.norm(X, 2) ** 2 / (len(y) * self.delta)
-
-    def get_lipschitz_sparse(self, X_data, X_indptr, X_indices, y):
-        """Compute coordinate-wise Lipschitz constants for sparse X."""
-        n_samples = len(y)
-        weight = max(self.quantile, 1 - self.quantile)
-        n_features = len(X_indptr) - 1
-        lipschitz = np.zeros(n_features, dtype=X_data.dtype)
-        for j in range(n_features):
-            nrm2 = 0.0
-            for idx in range(X_indptr[j], X_indptr[j + 1]):
-                nrm2 += X_data[idx] ** 2
-            lipschitz[j] = weight * nrm2 / (n_samples * self.delta)
-        return lipschitz
-
-    def get_global_lipschitz_sparse(self, X_data, X_indptr, X_indices, y):
-        """Compute global Lipschitz constant for sparse X."""
+    def gradient_scalar(self, X, y, w, Xw, j):
+        """Compute gradient w.r.t. w_j - following parent class pattern."""
         n_samples = len(y)
-        weight = max(self.quantile, 1 - self.quantile)
-        return weight * spectral_norm(
-            X_data, X_indptr, X_indices, n_samples
-        ) ** 2 / (n_samples * self.delta)
+        grad_j = 0.0
+        for i in range(n_samples):
+            residual = y[i] - Xw[i]
+            grad_j += -X[i, j] * self._grad_scalar(residual)
+        return grad_j / n_samples
+
+    def _grad_scalar(self, residual):
+        """Calculate gradient for a single residual."""
+        tau = self.quantile
+        delta = self.delta
+        r = residual
+
+        if r >= delta:
+            # Upper linear tail: r >= delta
+            return tau
+        elif r >= 0:
+            # Upper quadratic: 0 <= r < delta
+            return tau * r / delta
+        elif r > -delta:
+            # Lower quadratic: -delta < r < 0
+            return (1 - tau) * r / delta
+        else:
+            # Lower linear tail: r <= -delta
+            return tau - 1
+
+
+class SimpleQuantileRegressor(BaseEstimator, RegressorMixin):
+    """Simple quantile regression without progressive smoothing."""
+
+    def __init__(self, quantile=0.5, alpha=0.1, delta=0.1, max_iter=1000, tol=1e-4):
+        self.quantile = quantile
+        self.alpha = alpha
+        self.delta = delta
+        self.max_iter = max_iter
+        self.tol = tol
+
+    def fit(self, X, y):
+        """Fit using FISTA with fixed delta."""
+        X, y = check_X_y(X, y)
+
+        datafit = QuantileHuber(quantile=self.quantile, delta=self.delta)
+        penalty = L1(alpha=self.alpha)
+        solver = FISTA(max_iter=self.max_iter, tol=self.tol)
+
+        est = GeneralizedLinearEstimator(
+            datafit=datafit,
+            penalty=penalty,
+            solver=solver
+        )
+
+        est.fit(X, y)
+        self.coef_ = est.coef_
+
+        return self
 
-    def intercept_update_step(self, y, Xw):
-        return -np.mean(self.raw_grad(y, Xw))
+    def predict(self, X):
+        """Predict using the fitted model."""
+        X = check_array(X)
+        return X @ self.coef_