BUG & FEAT Fix convergence issues with Pinball loss on large datasets (Issue #276)

alpha_search.py

@@ -0,0 +1,36 @@
+import numpy as np
+from numpy.linalg import norm
+from skglm import Lasso
+from skglm.experimental import SqrtLasso
+
+np.random.seed(0)
+X = np.random.randn(10, 20)
+y = np.random.randn(10)
+y += 1
+
+n = len(y)
+
+alpha_max = norm(X.T @ y, ord=np.inf) / n
+
+alpha = alpha_max / 10
+
+lass = Lasso(alpha=alpha, fit_intercept=True, tol=1e-8).fit(X, y)
+w_lass = lass.coef_
+assert norm(w_lass) > 0
+
+scal = n / norm(y - lass.predict(X))
+
+sqrt = SqrtLasso(alpha=alpha * scal, fit_intercept=True, tol=1e-8).fit(X, y)
+
+print(norm(w_lass - sqrt.coef_))
+
+
+# diffs = []
+# alphas = np.linspace(16.07, 16.08, num=50)
+# for scal in alphas:
+#     sqrt = SqrtLasso(alpha=alpha * scal).fit(X, y)
+#     diffs.append(norm(w_lass - sqrt.coef_))
+
+# best = np.argmin(diffs)
+# print(alphas[best])
+# print(diffs[best])

doc/api.rst

@@ -107,3 +107,5 @@ Experimental
    Pinball
    SqrtQuadratic
    SqrtLasso
+   QuantileHuber
+   SmoothQuantileRegressor

doc/changes/0.5.rst

@@ -3,3 +3,4 @@
 Version 0.5 (in progress)
 -------------------------
 - Add support for fitting an intercept in :ref:`SqrtLasso <skglm.experimental.sqrt_lasso.SqrtLasso>` (PR: :gh:`298`)
+- Add :ref:`SmoothQuantileRegressor <skglm.experimental.smooth_quantile_regressor.SmoothQuantileRegressor>` for fast quantile regression with progressive smoothing (PR: :gh:`276`)

example.py

@@ -0,0 +1,52 @@
+import time
+import numpy as np
+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 sklearn.model_selection import train_test_split
+
+from numpy.linalg import norm
+
+
+def pinball_loss(y_true, y_pred, tau=0.5):
+    """Compute Pinball (quantile) loss."""
+    residuals = y_true - y_pred
+    return np.mean(np.where(residuals >= 0,
+                            tau * residuals,
+                            (1 - tau) * -residuals))
+
+
+# Test different problem sizes
+n_samples, n_features = 100, 100
+X, y = make_regression(n_samples=n_samples, n_features=n_features,
+                       noise=0.1, random_state=0)
+alpha = 0.01
+
+# Test different noise distributions
+
+# Quantiles to test
+tau = 0.3
+
+# Store results
+results = []
+
+X_train, X_test, y_train, y_test = train_test_split(
+    X, y, test_size=0.2, random_state=42
+)
+
+# scikit-learn QuantileRegressor
+qr = QuantileRegressor(quantile=tau, alpha=alpha, solver="highs")
+t0 = time.time()
+qr.fit(X_train, y_train)
+qr_time = time.time() - t0
+
+
+ours = SmoothQuantileRegressor(quantile=tau, alpha=alpha)
+t0 = time.time()
+ours.fit(X_train, y_train)
+ours_time = time.time() - t0
+
+
+print(ours.coef_ - qr.coef_)
+print(norm(ours.coef_ - qr.coef_) / norm(qr.coef_))