Issue with Lasso for zero weights

[idc9]

There is a numerical issue (i.e. the wrong answer is output) with celer.homotopy.celer_path when the user supplies a zero L1 weight vector. This use case comes up, for example, in reweighted Lasso problems (e.g. #185).

Relatedly when alpha=0 is specified there are warnings about convergence, however, the solution appears to be correct.

Reproducible experiment

I'm using: celer version 0.6, python version 3.6.12, and OSX Darwin.

from sklearn.datasets import make_regression
from sklearn.linear_model import Lasso, LinearRegression
import numpy as np

from celer.homotopy import celer_path

# make toy regression dataset
X, y = make_regression(n_samples=100, n_features=20, noise=1, random_state=324)

First let's do a sanity check to make sure celer and sklearn Lasso agree when alpha=1

est_lasso = Lasso(alpha=1, fit_intercept=False).fit(X, y)
coef_celer = celer_path(X=X, y=y, pb='lasso', alphas=[1])[1].reshape(-1)
print(np.linalg.norm(coef_celer - est_lasso.coef_))  

0.0008  # awesome they agree!

Bug with zero weight vector

est_ls = LinearRegression(fit_intercept=False).fit(X, y) # use sklearn to get baseline least square solution

# set the L1 weights to be the zero vector
weights = np.zeros(X.shape[1])

coef_celer = celer_path(X=X, y=y, pb='lasso',
                        alphas=[1], 
                        weights=weights)[1].reshape(-1)

I get the following warning

/Users/..../lib/python3.6/site-packages/celer/homotopy.py:290: RuntimeWarning: divide by zero encountered in true_divide
  theta /= scal

print(np.linalg.norm(coef_celer - est_ls.coef_)) 

160.6176631165807 # celer gives the wrong answer

Here celer gives the wrong answer!

Minor issue with alpha = 0

Next try to get the unpenalized least squares solution with celer (i.e. set alpha=0)

coef_celer = celer_path(X=X, y=y, pb='lasso', alphas=[0])[1].reshape(-1)

I get the following warning

!!! Inner solver did not converge at epoch 49999, gap: 7.51e+02 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06
!!! Inner solver did not converge at epoch 49999, gap: 4.73e-01 > 1.00e-06

print(np.linalg.norm(coef_celer - est_ls.coef_)) # good
4.172701492321924e-13 # awesome they agree!

celer gives the correct answer, but prints out convergence warnings. I don't know if this is a problem (e.g. I haven't done time tests to see if it's slowing down).

[mathurinm]

Hi Iain,
Thanks for the neat issue. In its current state, Celer relies heavily on the dual problem, and the dual variable theta must satisfy the constraint:
norm(X.T @ theta / weights, ord=np.inf) <= alpha. As you observed, this is an issue if some weights are 0, or if alpha is 0.

• in a push 2 days ago, I raised a ValueError is weights are not all > 0. weights equal to 0 used to encode that the feature should be ignored, this was bad for the user and has been changed to weight = np.inf
• celer is built for the Lasso (exploiting sparsity), so supporting alpha=0 is not on the menu for now. sklearn raises a warning in that case and I thought that we would, too.

Both of these can be addressed using a modification of celer implementing extrapolation in the primal as in https://github.com/mathurinm/andersoncd, which we'll do in the coming months

[mathurinm]

We have support for vanishing weights it in mathurinm/andersoncd#21 (in numba). Feedback is welcome.
I have to plan a bit but the goal is to integrate it to celer at some point. It requires some work to guarantee performance does not decrease, so it'll take a bit of time.
thanks for suggesting this features.