MNT - remove normalization by sqrt(n_samples) in Square root Lasso (#130)

Author: Badr-MOUFAD

skglm/experimental/sqrt_lasso.py

@@ -12,10 +12,11 @@
 
 
 class SqrtQuadratic(BaseDatafit):
-    """Square root quadratic datafit.
+    """Unnormalized square root quadratic datafit.
 
     The datafit reads::
-        ||y - Xw||_2 / sqrt(n_samples)
+
+        ||y - Xw||_2
     """
 
     def __init__(self):
@@ -29,7 +30,7 @@ def params_to_dict(self):
         return dict()
 
     def value(self, y, w, Xw):
-        return np.linalg.norm(y - Xw) / np.sqrt(len(y))
+        return np.linalg.norm(y - Xw)
 
     def raw_grad(self, y, Xw):
         """Compute gradient of datafit w.r.t ``Xw``.
@@ -45,12 +46,12 @@ def raw_grad(self, y, Xw):
         if norm_residuals < 1e-2 * norm(y):
             raise ValueError("SmallResidualException")
 
-        return minus_residual / (norm_residuals * np.sqrt(len(y)))
+        return minus_residual / norm_residuals
 
     def raw_hessian(self, y, Xw):
         """Diagonal matrix upper bounding the Hessian."""
         n_samples = len(y)
-        fill_value = 1 / (np.sqrt(n_samples) * norm(y - Xw))
+        fill_value = 1 / norm(y - Xw)
         return np.full(n_samples, fill_value)
 
 
@@ -59,7 +60,7 @@ class SqrtLasso(LinearModel, RegressorMixin):
 
     The optimization objective for square root Lasso is::
 
-        |y - X w||_2 / sqrt(n_samples) + alpha * ||w||_1
+        |y - X w||_2 + alpha * ||w||_1
 
     Parameters
     ----------
@@ -205,7 +206,8 @@ def _chambolle_pock_sqrt(X, y, alpha, max_iter=1000, obj_freq=10, verbose=False)
     """Apply Chambolle-Pock algorithm to solve square-root Lasso.
 
     The objective function is:
-        min_w ||Xw - y||_2/sqrt(n_samples) + alpha * ||w||_1.
+
+        min_w ||Xw - y||_2 + alpha * ||w||_1.
     """
     n_samples, n_features = X.shape
     # dual variable is z, primal is w
@@ -221,12 +223,12 @@ def _chambolle_pock_sqrt(X, y, alpha, max_iter=1000, obj_freq=10, verbose=False)
     sigma = 0.99 / L
 
     for t in range(max_iter):
-        w = ST_vec(w - tau * X.T @ (2 * z - z_old), alpha * np.sqrt(n_samples) * tau)
+        w = ST_vec(w - tau * X.T @ (2 * z - z_old), alpha * tau)
         z_old = z.copy()
         z[:] = proj_L2ball(z + sigma * (X @ w - y))
 
         if t % obj_freq == 0:
-            objs.append(norm(X @ w - y) / np.sqrt(n_samples) + alpha * norm(w, ord=1))
+            objs.append(norm(X @ w - y) + alpha * norm(w, ord=1))
             if verbose:
                 print(f"Iter {t}, obj {objs[-1]: .10f}")
 

skglm/experimental/tests/test_sqrt_lasso.py

@@ -9,7 +9,7 @@
 def test_alpha_max():
     n_samples, n_features = 50, 10
     X, y, _ = make_correlated_data(n_samples, n_features, random_state=0)
-    alpha_max = norm(X.T @ y, ord=np.inf) / (np.sqrt(n_samples) * norm(y))
+    alpha_max = norm(X.T @ y, ord=np.inf) / norm(y)
 
     sqrt_lasso = SqrtLasso(alpha=alpha_max).fit(X, y)
 
@@ -24,7 +24,7 @@ def test_vs_statsmodels():
     n_samples, n_features = 50, 10
     X, y, _ = make_correlated_data(n_samples, n_features, random_state=0)
 
-    alpha_max = norm(X.T @ y, ord=np.inf) / (np.sqrt(n_samples) * norm(y))
+    alpha_max = norm(X.T @ y, ord=np.inf) / norm(y)
     n_alphas = 3
     alphas = alpha_max * np.geomspace(1, 1e-2, n_alphas+1)[1:]
 
@@ -36,9 +36,10 @@ def test_vs_statsmodels():
     # fit statsmodels on path
     for i in range(n_alphas):
         alpha = alphas[i]
+        # statsmodels solves: ||y - Xw||_2 + alpha * ||w||_1 / sqrt(n_samples)
         model = linear_model.OLS(y, X)
         model = model.fit_regularized(method='sqrt_lasso', L1_wt=1.,
-                                      alpha=n_samples * alpha)
+                                      alpha=np.sqrt(n_samples) * alpha)
         coefs_statsmodels[i] = model.params
 
     np.testing.assert_almost_equal(coefs_skglm, coefs_statsmodels, decimal=4)
@@ -48,7 +49,7 @@ def test_prox_newton_cp():
     n_samples, n_features = 50, 10
     X, y, _ = make_correlated_data(n_samples, n_features, random_state=0)
 
-    alpha_max = norm(X.T @ y, ord=np.inf) / (np.sqrt(n_samples) * norm(y))
+    alpha_max = norm(X.T @ y, ord=np.inf) / norm(y)
     alpha = alpha_max / 10
     clf = SqrtLasso(alpha=alpha, tol=1e-12).fit(X, y)
     w, _, _ = _chambolle_pock_sqrt(X, y, alpha, max_iter=1000)