Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correctly handling the case λmax = 0. #53

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Correctly handling the case λmax = 0. #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3ccd5f0
Correctly handling the case λmax = 0.
barankarakus Sep 6, 2020
f4a3923
Update src/Lasso.jl
barankarakus Sep 20, 2020
5636fed
Replacing equality with approximate equality.
Sep 20, 2020
5ac04f7
Added test for case: zero variation in y.
Sep 20, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions src/Lasso.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,10 @@ const MAX_DEV_FRAC = 0.999
# Compute automatic λ values based on λmax and λminratio
function computeλ(λmax, λminratio, α, nλ)
λmax /= α
if isapprox(λmax, 0; atol=1e-10) # then assuming λmax = 0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is tricky because I think lambda is not unitless, so if it is small or not depends on the data given.
How does glmnet in R (or GLMNet.jl) handle this case?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reason I've changed the equality check to an isapprox() check is due to floating point arithmetic leading to a lambdamax that should actually be zero being very close to zero but non-zero instead. Simple example when this happens is a design matrix X with entries sampled from U[0, 1] and y a non-zero vector with identical entries.

I agree, the data could be such that lambdamax is genuinely very small but non-zero.

That said, I think it would be very rare to encounter such data in practice... especially since lambdamax (for the linear model) scales linearly with X and y, and we tend to standardise these.

I see two approaches going forward:

  1. Keep this check as is - the case where it would fail to produce correct output basically never occurs anyway.

  2. Revert back to the equality check. The real case in which the package failed was the case where lambdamax was exactly zero, anyway. Moreover, even if lambdamax should be zero but instead is a very small number, there is no major problem: the solver works very fast and it is clear from the output that every value of lambda yields zero active coefficients.

I'll leave it to you to decide 😃.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Additionally: I'm not sure how glmnet in R or Julia handles this.

@info "The penalized coefficients equal zero for all values of the regularization parameter λ."
return [0]
end
logλmax = log(λmax)
exp.(range(logλmax, stop=logλmax + log(λminratio), length=nλ))
end
Expand Down
40 changes: 22 additions & 18 deletions src/coordinate_descent.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -685,7 +685,7 @@ function StatsBase.fit!(path::RegularizationPath{S,T}; verbose::Bool=false, irls
niter = 0
if nλ == 0
i = 0
else
elseif i <= nλ # need this check because it is possible that autoλ is true and nλ is 1
while true # outer loop
obj = convert(T, Inf)
last_dev_ratio = dev_ratio
Expand Down Expand Up @@ -776,6 +776,7 @@ function StatsBase.fit!(path::RegularizationPath{S,T}; verbose::Bool=false, irls
end
end

i = min(i, nλ)
path.λ = path.λ[1:i]
path.pct_dev = pct_dev[1:i]
path.coefs = coefs[:, 1:i]
Expand Down Expand Up @@ -819,29 +820,32 @@ function StatsBase.fit!(path::RegularizationPath{S,T}; verbose::Bool=false,
i = 1
end

while true # outer loop
last_dev_ratio = dev_ratio
curλ = λ[i]
if i <= nλ # need this check because it is possible that autoλ is true and nλ is 1
while true # outer loop
last_dev_ratio = dev_ratio
curλ = λ[i]

# Run coordinate descent
niter += cdfit!(newcoef, cd, curλ, criterion)
# Run coordinate descent
niter += cdfit!(newcoef, cd, curλ, criterion)

dev_ratio = cd.dev/nulldev
pct_dev[i] = 1 - dev_ratio
addcoefs!(coefs, newcoef, i)
b0s[i] = intercept(newcoef, cd)
dev_ratio = cd.dev/nulldev
pct_dev[i] = 1 - dev_ratio
addcoefs!(coefs, newcoef, i)
b0s[i] = intercept(newcoef, cd)

# Test whether we should continue
if i == nλ || (stopearly && autoλ && (last_dev_ratio - dev_ratio < MIN_DEV_FRAC_DIFF ||
pct_dev[i] > MAX_DEV_FRAC))
break
end
# Test whether we should continue
if i == nλ || (stopearly && autoλ && (last_dev_ratio - dev_ratio < MIN_DEV_FRAC_DIFF ||
pct_dev[i] > MAX_DEV_FRAC))
break
end

verbose && println("$i: λ=$curλ, pct_dev=$(pct_dev[i])")
poststep(path, cd, i, newcoef)
i += 1
verbose && println("$i: λ=$curλ, pct_dev=$(pct_dev[i])")
poststep(path, cd, i, newcoef)
i += 1
end
end

i = min(i, nλ)
path.λ = path.λ[1:i]
path.pct_dev = pct_dev[1:i]
path.coefs = coefs[:, 1:i]
Expand Down
18 changes: 18 additions & 0 deletions test/lasso.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,6 +154,24 @@ end
end
end

# Test case with zero variation in y is handled correctly
function zero_variation_test()
X = [
0.5472502169628388 0.37660447632078875 0.06669114126498532 0.4950818154768257;
0.5142931961160688 0.520205941129849 0.4052730635141131 0.6700530909562794;
0.5831846867316071 0.3174143498124731 0.772131243876973 0.03386847158881201;
0.8802489459954292 0.6742158685234003 0.3849775799923969 0.7773264968613842;
0.9216786846192617 0.7888303438159934 0.09788865152005011 0.34950775139369905
]
y = 0.2937233091452627 .+ zeros(size(X, 1))
path = fit(LassoPath, X, y)
(path.λ == eltype(path.λ)[0]) || return false
(length(path.coefs.nzval) == 0) || return false
return true
end

@test zero_variation_test() == true
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe use @test_log instead?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, any idea why the tests stopped working in julia v1.0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe use @test_log instead?

I agree. Will implement tomorrow.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, any idea why the tests stopped working in julia v1.0?

Unfortunately nope!


# Test for sparse matrices

# @testset "LassoPath Zero in" begin
Expand Down