Formula for dual gap of elastic net in coordinate descent solver

### Describe the bug
The computation of the dual gap for the elastic net in the coordinate descent solver (`enet_coordinate_descent`) might be wrong.

The elastic net minimizes
```
Primal(w) = (1/2) * ||y - X w||_2^2 + alpha * ||w||_1 + beta/2 * ||w||_2^2
```

### Lasso
For the pure Lasso, i.e. `beta=0`, the dual becomes [1]
```
Dual(nu) = -1/2 * ||nu||_2^2 - nu'y   if ||X'nu||_∞ = max(abs(X'nu)) <= alpha   else -∞
```
with `nu = Xw - y` (=minus the residual `R`) as possible dual point. This yields the Lasso dual gap
```
Primal(w) - Dual(nu)
```
In case of `||X'nu||_∞ > alpha`, one uses a down-scaled variable `nu` in the dual.

### Elastic Net
For the elastic net, the dual, see section 5.2.3 in [2], becomes
```
Dual(nu) = -1/2 * ||nu||_2^2 - nu'y - 1/(2*beta) * sum_j (|X_j'nu| - alpha)_+^2
```
with `X_j` the j-th column of `X` and `x_+ = max(0, x)` the positive part.

### References
[1] Kim, Seung-Jean et al. “An Interior-Point Method for Large-Scale $\ell_1$-Regularized Least Squares.” IEEE Journal of Selected Topics in Signal Processing 1 (2007): 606-617. [pdf link](http://stanford.edu/~boyd/papers/pdf/l1_ls.pdf)
[2] Mendler-Dünner, Celestine et al. “Primal-Dual Rates and Certificates.” ICML (2016). [arxiv link](https://arxiv.org/pdf/1602.05205.pdf)

### Expected Results

```python
if beta > =0:
    R = y - X @ w  # note: R = -nu
    R_norm2 = R @ R
    gap = R_norm2 - R @ y + alpha * l1_norm + beta/2 * l2_norm
    gap += 1/(2*beta) * np.sum(np.max(0, np.abs(X.T @ R) - alpha)**2)
```
At least a test for it, similar to the pure Lasso case in https://github.com/scikit-learn/scikit-learn/blob/2c2f31d68c21b7647557b2f776e86c05954d80bf/sklearn/linear_model/tests/test_coordinate_descent.py#L292

### Actual Results

I start to be rattled by the `-beta * w` term in line 205.
https://github.com/scikit-learn/scikit-learn/blob/2c2f31d68c21b7647557b2f776e86c05954d80bf/sklearn/linear_model/_cd_fast.pyx#L205-L236

COMMENT:
That's a great code snippet!

- `n_samples=500`, `n_features=2000`, `alpha = alpha_max / 30`, `beta = 0.1 * alpha`
  1583 zero coefficients `np.sum(w == 0)`
  ![image](https://user-images.githubusercontent.com/15324633/158613117-99387dce-a4d6-4209-853d-fe521415a918.png)
- `n_samples=500`, `n_features=2000`, `alpha = alpha_max / 30`, **`beta = 1.0 * alpha`**
   ![image](https://user-images.githubusercontent.com/15324633/158614549-87f7c700-ac96-4a67-b714-91cb24aa2b2e.png)
- `n_samples=500`, `n_features=2000`, **`alpha = alpha_max / 300`**, `beta = 0.1 * alpha`
  ![image](https://user-images.githubusercontent.com/15324633/158613831-f4cad3d9-2f50-4544-85a1-457706ff3f2c.png)
- `n_samples=500`, `n_features=2000`, **`alpha = alpha_max / 300`**, **`beta = 1.0 * alpha`**
  ![image](https://user-images.githubusercontent.com/15324633/158614704-18602f42-2db6-4735-bb44-386c742b7330.png)
- **`n_samples=5000`**, **`n_features=200`**, `alpha = alpha_max / 30`, **`beta = 0.1 * alpha`
  ![image](https://user-images.githubusercontent.com/15324633/158615055-7860c105-3865-4de6-b819-4d84d6ea077e.png)


**Summary:** It depends.
At the first few iterations, the Lasso formulation always wins. For later iterations, the unconstraint formulation seems favorably in most settings.

**Proposition:** Could it make sense to take the minimum of both versions? Most ingredients are already computed anyway.

Labels: Bug, module:linear_model