Repository: glm-tools/pyglmnet

Issue:
failing comparison of elastic net with scikit-learn
I tried to compare pyglmnet.GLM with ElasticNet from scikit-learn and could not get it work. Code to reproduce:

```python
import numpy as np
from sklearn.datasets.samples_generator import make_regression
from sklearn.linear_model import ElasticNet
from pyglmnet import GLM

def rmse(a, b):
    return np.sqrt(np.mean((a - b) ** 2))

X, Y, coef_ = make_regression(
    n_samples=1000, n_features=1000,
    noise=0.1, n_informative=10, coef=True,
    random_state=42)

alpha = 0.1
l1_ratio=0.5

sk = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, tol = 1e-5).fit(X, Y)
pg = GLM(distr='gaussian', alpha=l1_ratio, reg_lambda=alpha, solver='cdfast', tol = 1e-5).fit(X, Y)
print('in-sample rmse sklearn = {}, pyglmnet = {}'.format(rmse(Y, sk.predict(X)), rmse(Y, pg.predict(X))))
```

Result:
in-sample rmse sklearn = 12.756054997216442, pyglmnet = 161.03496460055504

Discussion:
peterfoley605: re. flake8 errors, I addressed everything that's touched by this PR (using `git diff main/master -U0  | flake8 --diff -`). The other flake8 issues are leftovers from `master`.
codecov-io: # [Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=h1) Report
> Merging [#278](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=desc) into [master](https://codecov.io/gh/glm-tools/pyglmnet/commit/09fdad1db0daed804a864ba6bdf36e3d20016e79?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/glm-tools/pyglmnet/pull/278/graphs/tree.svg?width=650&token=kkw3ydHRLW&height=150&src=pr)](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #278   +/-   ##
=======================================
  Coverage   73.74%   73.74%           
=======================================
  Files           4        4           
  Lines         617      617           
  Branches      124      124           
=======================================
  Hits          455      455           
- Misses        123      124    +1     
+ Partials       39       38    -1
```


| [Impacted Files](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pyglmnet/pyglmnet.py](https://codecov.io/gh/glm-tools/pyglmnet/pull/278/diff?src=pr&el=tree#diff-cHlnbG1uZXQvcHlnbG1uZXQucHk=) | `79.67% <100%> (ø)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=footer). Last update [09fdad1...5175d24](https://codecov.io/gh/glm-tools/pyglmnet/pull/278?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

jasmainak: @peterfoley do you mind if I push a commit to conform to the coding style so we can merge this?
peterfoley605: @jasmainak I just enabled "allow edits from maintainers" so you can push commits.
peterfoley605: @jasmainak I massively pared down this PR since the changes other than cdfast caching were resolved by other PRs over the past few months. Could you re-review?

As a reminder for why `z` from `_cdfast` shouldn't be cached outside `_cdfast` itself -- immediately after `_cdfast` is called by `fit()`, `beta` is updated by `self._prox`, which makes the cached `z` inconsistent with `beta`.

Patch:
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
index fe47a2d..565e172 100644
--- a/doc/whats_new.rst
+++ b/doc/whats_new.rst
@@ -23,6 +23,7 @@ Changelog
 BUG
 ~~~
 
+    - Fixed z cache inconsistency in cdfast when alpha and reg_lambda are nonzero by `Peter Foley`_. 
     - Fixed incorrect usage of the random_state parameter by `Giovanni De Toni`_.
     - Fixed incorrect proximal operator for group lasso by `Yu Umegaki`_.
     - Changed stopping criteria for convergence (a threshold on the change in objective value) which
@@ -71,3 +72,4 @@ Changelog
 .. _Giovanni De Toni: https://github.com/geektoni
 .. _Beibin Li: https://github.com/BeibinLi
 .. _Scott Otterson: https://github.com/notuntoward
+.. _Peter Foley: https://github.com/peterfoley
diff --git a/pyglmnet/pyglmnet.py b/pyglmnet/pyglmnet.py
index eb53fb4..e34eac2 100644
--- a/pyglmnet/pyglmnet.py
+++ b/pyglmnet/pyglmnet.py
@@ -615,7 +615,7 @@ class GLM(BaseEstimator):
 
             return result
 
-    def _cdfast(self, X, y, z, ActiveSet, beta, rl):
+    def _cdfast(self, X, y, ActiveSet, beta, rl):
         """
         Perform one cycle of Newton updates for all coordinates.
 
@@ -627,9 +627,6 @@ class GLM(BaseEstimator):
         y : array
             Labels to the data
             n_samples x 1
-        z:  array
-            n_samples x 1
-            beta[0] + X * beta[1:]
         ActiveSet: array
             n_features + 1 x 1
             Active set storing which betas are non-zero
@@ -644,13 +641,11 @@ class GLM(BaseEstimator):
         beta: array
             (n_features + 1) x 1
             Updated parameters
-        z: array
-            beta[0] + X * beta[1:]
-            (n_features + 1) x 1
         """
         n_samples = X.shape[0]
         n_features = X.shape[1]
         reg_scale = rl * (1 - self.alpha)
+        z = beta[0] + np.dot(X, beta[1:])
 
         for k in range(0, n_features + 1):
             # Only update parameters in active set
@@ -677,7 +672,7 @@ class GLM(BaseEstimator):
                 # Update parameters, z
                 update = 1. / hk * gk
                 beta[k], z = beta[k] - update, z - update * xk
-        return beta, z
+        return beta
 
     def fit(self, X, y):
         """The fit function.
@@ -743,7 +738,6 @@ class GLM(BaseEstimator):
 
         if self.solver == 'cdfast':
             ActiveSet = np.ones(n_features + 1)     # init active set
-            z = beta[0] + np.dot(X, beta[1:])       # cache z
 
         # Iterative updates
         for t in range(0, self.max_iter):
@@ -761,8 +755,8 @@ class GLM(BaseEstimator):
 
             elif self.solver == 'cdfast':
                 beta_old = deepcopy(beta)
-                beta, z = \
-                    self._cdfast(X, y, z, ActiveSet, beta, reg_lambda)
+                beta = \
+                    self._cdfast(X, y, ActiveSet, beta, reg_lambda)
                 # Converged if the norm(update) < tol
                 if (t > 1) and (np.linalg.norm(beta - beta_old) < tol):
                     msg = ('\tConverged in {0:d} iterations'.format(t))


Test patch:
diff --git a/tests/test_pyglmnet.py b/tests/test_pyglmnet.py
index 042ae00..a6536d8 100644
--- a/tests/test_pyglmnet.py
+++ b/tests/test_pyglmnet.py
@@ -10,6 +10,7 @@ from scipy.optimize import approx_fprime
 from sklearn.datasets import make_regression
 from sklearn.preprocessing import StandardScaler
 from sklearn.model_selection import GridSearchCV, cross_val_score, KFold
+from sklearn.linear_model import ElasticNet
 
 from pyglmnet import (GLM, GLMCV, _grad_L2loss, _L2loss, simulate_glm,
                       _gradhess_logloss_1d, _loss, datasets, ALLOWED_DISTRS)
@@ -150,7 +151,8 @@ def test_group_lasso():
 
 
 @pytest.mark.parametrize("distr", ALLOWED_DISTRS)
-def test_glmnet(distr):
+@pytest.mark.parametrize("reg_lambda", [0.0, 0.1])
+def test_glmnet(distr, reg_lambda):
     """Test glmnet."""
     raises(ValueError, GLM, distr='blah')
     raises(ValueError, GLM, distr='gaussian', max_iter=1.8)
@@ -182,14 +184,13 @@ def test_glmnet(distr):
                                sample=False)
 
         alpha = 0.
-        reg_lambda = 0.
         loss_trace = list()
+        eta = 2.0
+        group = None
+        Tau = None
 
         def callback(beta):
             Tau = None
-            eta = 2.0
-            group = None
-
             loss_trace.append(
                 _loss(distr, alpha, Tau, reg_lambda,
                       X_train, y_train, eta, group, beta))
@@ -205,13 +206,15 @@ def test_glmnet(distr):
         # verify loss decreases
         assert(np.all(np.diff(loss_trace) <= 1e-7))
 
-        # verify loss at convergence = loss when beta=beta_
-        l_true = _loss(distr, 0., np.eye(beta.shape[0]), 0.,
-                       X_train, y_train, 2.0, None,
-                       np.concatenate(([beta0], beta)))
-        assert_allclose(loss_trace[-1], l_true, rtol=1e-4, atol=1e-5)
-        # beta=beta_ when reg_lambda = 0.
-        assert_allclose(beta, glm.beta_, rtol=0.05, atol=1e-2)
+        # true loss and beta should be recovered when reg_lambda == 0
+        if reg_lambda == 0.:
+            # verify loss at convergence = loss when beta=beta_
+            l_true = _loss(distr, alpha, Tau, reg_lambda,
+                           X_train, y_train, eta, group,
+                           np.concatenate(([beta0], beta)))
+            assert_allclose(loss_trace[-1], l_true, rtol=1e-4, atol=1e-5)
+            # beta=beta_ when reg_lambda = 0.
+            assert_allclose(beta, glm.beta_, rtol=0.05, atol=1e-2)
         betas_.append(glm.beta_)
 
         y_pred = glm.predict(X_train)
@@ -293,6 +296,32 @@ def test_cv():
     glmcv.fit(X, y)
 
 
+@pytest.mark.parametrize("solver", ['batch-gradient', 'cdfast'])
+def test_compare_sklearn(solver):
+    """Test results against sklearn."""
+
+    def rmse(a, b):
+        return np.sqrt(np.mean((a - b) ** 2))
+
+    X, Y, coef_ = make_regression(
+        n_samples=1000, n_features=1000,
+        noise=0.1, n_informative=10, coef=True,
+        random_state=42)
+
+    alpha = 0.1
+    l1_ratio = 0.5
+
+    clf = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, tol=1e-5)
+    clf.fit(X, Y)
+    glm = GLM(distr='gaussian', alpha=l1_ratio, reg_lambda=alpha,
+              solver=solver, tol=1e-5, max_iter=70)
+    glm.fit(X, Y)
+
+    y_sk = clf.predict(X)
+    y_pg = glm.predict(X)
+    assert abs(rmse(Y, y_sk) - rmse(Y, y_pg)) < 1.0
+
+
 @pytest.mark.parametrize("distr", ALLOWED_DISTRS)
 def test_cdfast(distr):
     """Test all functionality related to fast coordinate descent."""
@@ -343,10 +372,8 @@ def test_cdfast(distr):
 
     # test cdfast
     ActiveSet = np.ones(n_features + 1)
-    beta_ret, z_ret = glm._cdfast(X, y, z,
-                                  ActiveSet, beta_, glm.reg_lambda)
+    beta_ret = glm._cdfast(X, y, ActiveSet, beta_, glm.reg_lambda)
     assert(beta_ret.shape == beta_.shape)
-    assert(z_ret.shape == z.shape)
 
 
 def test_fetch_datasets():