From 3232c370b4ff8a9c6f607f00c829ff9002a9c287 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Tue, 30 Nov 2021 16:17:03 +0100 Subject: [PATCH 01/34] Fix bug in glmnet solver --- solvers/glmnet.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 7431a67c..17b95a65 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -42,14 +42,15 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] def run(self, n_iter): - fit_dict = {"lambda.min.ratio": self.lmbd / self.lmbd_max} + fit_dict = {"lambda": self.lmbd} + glmnet_fit = self.glmnet(self.X, self.y, intercept=False, standardize=False, maxit=n_iter, thresh=1e-14, **fit_dict) results = dict(zip(glmnet_fit.names, list(glmnet_fit))) as_matrix = robjects.r['as'] coefs = np.array(as_matrix(results["beta"], "matrix")) - self.w = coefs[:, -1] + self.w = coefs.flatten() def get_result(self): return self.w From df67ce628b68db36e321a8c626ade614121b583a Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 1 Dec 2021 11:06:54 +0100 Subject: [PATCH 02/34] add debug script --- test_glmnet.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 test_glmnet.py diff --git a/test_glmnet.py b/test_glmnet.py new file mode 100644 index 00000000..6878c354 --- /dev/null +++ b/test_glmnet.py @@ -0,0 +1,36 @@ +import numpy as np +from rpy2 import robjects +from rpy2.robjects import numpy2ri +from rpy2.robjects.packages import importr +import celer +from numpy.linalg import norm + +numpy2ri.activate() +glmnet = importr('glmnet') +as_matrix = robjects.r['as'] + +n = 10 +p = 30 + +np.random.seed(1532) + +y = np.random.randn(n) +X = np.random.randn(n, p) + +lmbd_max = np.max(np.abs(X.T @ y)) / n + +for ratio in [0.1, 1e-3]: + lmbd = ratio * lmbd_max + + fit_dict2 = {"lambda": lmbd} + fit2 = glmnet.glmnet(X, y, intercept=False, + standardize=False, **fit_dict2) + results2 = dict(zip(fit2.names, list(fit2))) + + beta = np.squeeze(np.array(as_matrix(results2["beta"], "matrix"))) + + clf = celer.Lasso(fit_intercept=False, alpha=lmbd, tol=1e-6).fit(X, y) + print(f"lambda = {ratio}*lambda_max") + print(beta) + print(clf.coef_) + print(f"norm of difference: {norm(beta - clf.coef_)}") From 6aa76de981576ca295d703e48257d63a08e561f0 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 1 Dec 2021 11:14:12 +0100 Subject: [PATCH 03/34] add computation of objective functions --- test_glmnet.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test_glmnet.py b/test_glmnet.py index 6878c354..5e4b191f 100644 --- a/test_glmnet.py +++ b/test_glmnet.py @@ -19,6 +19,12 @@ lmbd_max = np.max(np.abs(X.T @ y)) / n + +def p_obj(X, y, w, lmbd): + R = y - X @ w + return np.mean(R ** 2) / 2. + lmbd * norm(w, ord=1) + + for ratio in [0.1, 1e-3]: lmbd = ratio * lmbd_max @@ -34,3 +40,5 @@ print(beta) print(clf.coef_) print(f"norm of difference: {norm(beta - clf.coef_)}") + print(f"celer obj: {p_obj(X, y, clf.coef_, lmbd)}") + print(f"GLMnet obj: {p_obj(X, y, beta, lmbd)}") From 47ed26502c8bf0fb2ddc7883b87dbbaee215d26e Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 1 Dec 2021 11:16:35 +0100 Subject: [PATCH 04/34] test that for lambda just below lambda max both solutions are non zero --- test_glmnet.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test_glmnet.py b/test_glmnet.py index 5e4b191f..efbcdc95 100644 --- a/test_glmnet.py +++ b/test_glmnet.py @@ -25,7 +25,8 @@ def p_obj(X, y, w, lmbd): return np.mean(R ** 2) / 2. + lmbd * norm(w, ord=1) -for ratio in [0.1, 1e-3]: +for ratio in [1, 0.95, 1e-3]: + print("#" * 80) lmbd = ratio * lmbd_max fit_dict2 = {"lambda": lmbd} From 4d2464802a5c10aaa04d563f9cd2cab536d71906 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 1 Dec 2021 11:26:54 +0100 Subject: [PATCH 05/34] Drop unnecessary computation of lambda_max --- solvers/glmnet.py | 1 - 1 file changed, 1 deletion(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 17b95a65..d790d7e0 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -38,7 +38,6 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept - self.lmbd_max = np.max(np.abs(X.T @ y)) self.glmnet = robjects.r['glmnet'] def run(self, n_iter): From fed8fa451c19873d5db0875444e48b1c7079fcb8 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 1 Dec 2021 11:27:36 +0100 Subject: [PATCH 06/34] Harmonize lambda scaling for glmnet solver --- solvers/glmnet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index d790d7e0..5cdb5326 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -41,7 +41,7 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] def run(self, n_iter): - fit_dict = {"lambda": self.lmbd} + fit_dict = {"lambda": self.lmbd / self.X.shape[0]} glmnet_fit = self.glmnet(self.X, self.y, intercept=False, standardize=False, maxit=n_iter, From 790607a13dd53ade500d0ce74f509e2bd7501899 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Mon, 6 Dec 2021 13:51:33 +0100 Subject: [PATCH 07/34] Switch benchopt method to use tolerance instead --- solvers/glmnet.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 5cdb5326..4d95fec4 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -24,7 +24,7 @@ class Solver(BaseSolver): 'for generalized linear models via coordinate descent", ' 'J. Stat. Softw., vol. 33, no. 1, pp. 1-22, NIH Public Access (2010)' ] - stop_strategy = 'iteration' + stop_strategy = 'tolerance' support_sparse = False def skip(self, X, y, lmbd, fit_intercept): @@ -40,12 +40,14 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] - def run(self, n_iter): - fit_dict = {"lambda": self.lmbd / self.X.shape[0]} + def run(self, tol): + fit_dict = {"lambda": self.lmbd / len(self.y)} + + print(tol) glmnet_fit = self.glmnet(self.X, self.y, intercept=False, - standardize=False, maxit=n_iter, - thresh=1e-14, **fit_dict) + standardize=False, maxit=1_000_000, + thresh=tol / len(self.y), **fit_dict) results = dict(zip(glmnet_fit.names, list(glmnet_fit))) as_matrix = robjects.r['as'] coefs = np.array(as_matrix(results["beta"], "matrix")) From 8f947d8710fc1360404bdd62a45dd836b7cc41f4 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Mon, 6 Dec 2021 15:38:12 +0100 Subject: [PATCH 08/34] tweaks to test_glmnet --- test_glmnet.py | 89 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 65 insertions(+), 24 deletions(-) diff --git a/test_glmnet.py b/test_glmnet.py index efbcdc95..366ef66d 100644 --- a/test_glmnet.py +++ b/test_glmnet.py @@ -1,45 +1,86 @@ +import sys +from timeit import default_timer as timer + +import benchopt +from benchopt.datasets import make_correlated_data import numpy as np +from numpy.linalg import norm from rpy2 import robjects from rpy2.robjects import numpy2ri from rpy2.robjects.packages import importr -import celer -from numpy.linalg import norm +import matplotlib.pyplot as plt +from celer import Lasso +import time +# Setup the system to allow rpy2 running numpy2ri.activate() glmnet = importr('glmnet') + as_matrix = robjects.r['as'] -n = 10 -p = 30 +n = 100 +p = 500 np.random.seed(1532) y = np.random.randn(n) -X = np.random.randn(n, p) +x, y, _ = make_correlated_data(n, p, random_state=0) + +lmbd_max = np.max(np.abs(x.T @ y)) + +lmbd = lmbd_max * 0.01 + +fit_dict = {"lambda": lmbd / n} + +tols = np.geomspace(1e-1, 1e-16, 11) + + +def compute_gap(x, y, beta, lmbd): + # compute residuals + diff = y - x.dot(beta) + + # compute primal objective and duality gap + p_obj = .5 * diff.dot(diff) + lmbd * np.abs(beta).sum() + theta = diff / lmbd + theta /= norm(x.T @ theta, ord=np.inf) + d_obj = (norm(y)**2 / 2. - lmbd**2 * norm(y / lmbd - theta)**2 / 2) + + return p_obj - d_obj + + +n_tols = len(tols) + +times = np.empty(n_tols) +gaps = np.empty(n_tols) -lmbd_max = np.max(np.abs(X.T @ y)) / n +for i, tol in enumerate(tols): + print(i, tol) + t0 = timer() + fit = glmnet.glmnet(x, + y, + intercept=False, + standardize=False, + **fit_dict, + thresh=tol, + maxit=10_000_000) -def p_obj(X, y, w, lmbd): - R = y - X @ w - return np.mean(R ** 2) / 2. + lmbd * norm(w, ord=1) + times[i] = timer() - t0 + results = dict(zip(fit.names, list(fit))) + coefs = np.array(as_matrix(results["beta"], "matrix")) + beta = coefs.flatten() -for ratio in [1, 0.95, 1e-3]: - print("#" * 80) - lmbd = ratio * lmbd_max + gaps[i] = compute_gap(x, y, beta, lmbd) - fit_dict2 = {"lambda": lmbd} - fit2 = glmnet.glmnet(X, y, intercept=False, - standardize=False, **fit_dict2) - results2 = dict(zip(fit2.names, list(fit2))) +plt.close('all') +plt.semilogy(times, gaps) +plt.xlabel("time (s)") +plt.ylabel("duality gap") +plt.show(block=False) - beta = np.squeeze(np.array(as_matrix(results2["beta"], "matrix"))) - clf = celer.Lasso(fit_intercept=False, alpha=lmbd, tol=1e-6).fit(X, y) - print(f"lambda = {ratio}*lambda_max") - print(beta) - print(clf.coef_) - print(f"norm of difference: {norm(beta - clf.coef_)}") - print(f"celer obj: {p_obj(X, y, clf.coef_, lmbd)}") - print(f"GLMnet obj: {p_obj(X, y, beta, lmbd)}") +t0 = time.time() +clf = Lasso(alpha=lmbd/n, fit_intercept=False, tol=1e-14, verbose=1).fit(x, y) +t1 = time.time() +print(t1 - t0) From 1be208a9b3cb7f9170bc3b6ece0f83ec19efeeef Mon Sep 17 00:00:00 2001 From: mathurinm Date: Mon, 6 Dec 2021 17:37:00 +0100 Subject: [PATCH 09/34] change glmnet criterion to use a higher patience --- solvers/glmnet.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 4d95fec4..9e1256d2 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -1,5 +1,6 @@ from benchopt import BaseSolver from benchopt import safe_import_context +from benchopt.stopping_criterion import SufficientProgressCriterion with safe_import_context() as import_ctx: @@ -24,9 +25,11 @@ class Solver(BaseSolver): 'for generalized linear models via coordinate descent", ' 'J. Stat. Softw., vol. 33, no. 1, pp. 1-22, NIH Public Access (2010)' ] - stop_strategy = 'tolerance' support_sparse = False + stopping_criterion = SufficientProgressCriterion( + patience=5, eps=1e-38, strategy='tolerance') + def skip(self, X, y, lmbd, fit_intercept): # XXX - glmnet support intercept, adapt the API if fit_intercept: @@ -43,8 +46,6 @@ def set_objective(self, X, y, lmbd, fit_intercept): def run(self, tol): fit_dict = {"lambda": self.lmbd / len(self.y)} - print(tol) - glmnet_fit = self.glmnet(self.X, self.y, intercept=False, standardize=False, maxit=1_000_000, thresh=tol / len(self.y), **fit_dict) From 46f6f0436c3da558405ed7184505f387d14d41b3 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Tue, 7 Dec 2021 14:01:53 +0100 Subject: [PATCH 10/34] faster tol decrease for glmnet, try to have common initial point --- solvers/glmnet.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 9e1256d2..990dc5a7 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -1,5 +1,5 @@ -from benchopt import BaseSolver -from benchopt import safe_import_context +from benchopt import BaseSolver, safe_import_context +from benchopt.runner import INFINITY from benchopt.stopping_criterion import SufficientProgressCriterion @@ -44,11 +44,12 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] def run(self, tol): + maxit = 0 if tol == INFINITY else 1_000_000 fit_dict = {"lambda": self.lmbd / len(self.y)} glmnet_fit = self.glmnet(self.X, self.y, intercept=False, - standardize=False, maxit=1_000_000, - thresh=tol / len(self.y), **fit_dict) + standardize=False, maxit=maxit, + thresh=tol ** 2.3, **fit_dict) results = dict(zip(glmnet_fit.names, list(glmnet_fit))) as_matrix = robjects.r['as'] coefs = np.array(as_matrix(results["beta"], "matrix")) From 68a379cab9666a811371c331b2bdb41dd0e35118 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Tue, 7 Dec 2021 14:40:50 +0100 Subject: [PATCH 11/34] increase patience again --- solvers/glmnet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 990dc5a7..95d427d8 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -28,7 +28,7 @@ class Solver(BaseSolver): support_sparse = False stopping_criterion = SufficientProgressCriterion( - patience=5, eps=1e-38, strategy='tolerance') + patience=7, eps=1e-38, strategy='tolerance') def skip(self, X, y, lmbd, fit_intercept): # XXX - glmnet support intercept, adapt the API From c57994db373fdf8f9d1a98963a79a3c50b96c681 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 8 Dec 2021 09:33:45 +0100 Subject: [PATCH 12/34] get common starting point, comment code --- solvers/glmnet.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 95d427d8..8b058446 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -44,12 +44,24 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] def run(self, tol): - maxit = 0 if tol == INFINITY else 1_000_000 + # even if maxit=0, glmnet can return non zero coefficients. To get the + # initial point on the curve, we set tol=0: this way, glmnet + # convergence check fails, and rather than returning the current + # iterate, it returns a 0 model. + if tol == INFINITY: + maxit = 0 + thresh = 0 + else: + maxit = 1_000_000 + # we need thresh to decay fast, otherwise the objective curve can + # plateau before convergence + thresh = tol ** 1.8 + fit_dict = {"lambda": self.lmbd / len(self.y)} glmnet_fit = self.glmnet(self.X, self.y, intercept=False, standardize=False, maxit=maxit, - thresh=tol ** 2.3, **fit_dict) + thresh=thresh, **fit_dict) results = dict(zip(glmnet_fit.names, list(glmnet_fit))) as_matrix = robjects.r['as'] coefs = np.array(as_matrix(results["beta"], "matrix")) From 43723de693f26cac7d244d44bc7a88b0aa72fd8d Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 8 Dec 2021 12:42:27 +0100 Subject: [PATCH 13/34] CLN remove test_glmnet.py --- test_glmnet.py | 86 -------------------------------------------------- 1 file changed, 86 deletions(-) delete mode 100644 test_glmnet.py diff --git a/test_glmnet.py b/test_glmnet.py deleted file mode 100644 index 366ef66d..00000000 --- a/test_glmnet.py +++ /dev/null @@ -1,86 +0,0 @@ -import sys -from timeit import default_timer as timer - -import benchopt -from benchopt.datasets import make_correlated_data -import numpy as np -from numpy.linalg import norm -from rpy2 import robjects -from rpy2.robjects import numpy2ri -from rpy2.robjects.packages import importr -import matplotlib.pyplot as plt -from celer import Lasso -import time - -# Setup the system to allow rpy2 running -numpy2ri.activate() -glmnet = importr('glmnet') - -as_matrix = robjects.r['as'] - -n = 100 -p = 500 - -np.random.seed(1532) - -y = np.random.randn(n) -x, y, _ = make_correlated_data(n, p, random_state=0) - -lmbd_max = np.max(np.abs(x.T @ y)) - -lmbd = lmbd_max * 0.01 - -fit_dict = {"lambda": lmbd / n} - -tols = np.geomspace(1e-1, 1e-16, 11) - - -def compute_gap(x, y, beta, lmbd): - # compute residuals - diff = y - x.dot(beta) - - # compute primal objective and duality gap - p_obj = .5 * diff.dot(diff) + lmbd * np.abs(beta).sum() - theta = diff / lmbd - theta /= norm(x.T @ theta, ord=np.inf) - d_obj = (norm(y)**2 / 2. - lmbd**2 * norm(y / lmbd - theta)**2 / 2) - - return p_obj - d_obj - - -n_tols = len(tols) - -times = np.empty(n_tols) -gaps = np.empty(n_tols) - -for i, tol in enumerate(tols): - print(i, tol) - t0 = timer() - - fit = glmnet.glmnet(x, - y, - intercept=False, - standardize=False, - **fit_dict, - thresh=tol, - maxit=10_000_000) - - times[i] = timer() - t0 - - results = dict(zip(fit.names, list(fit))) - coefs = np.array(as_matrix(results["beta"], "matrix")) - beta = coefs.flatten() - - gaps[i] = compute_gap(x, y, beta, lmbd) - -plt.close('all') -plt.semilogy(times, gaps) -plt.xlabel("time (s)") -plt.ylabel("duality gap") -plt.show(block=False) - - -t0 = time.time() -clf = Lasso(alpha=lmbd/n, fit_intercept=False, tol=1e-14, verbose=1).fit(x, y) -t1 = time.time() -print(t1 - t0) From 5fa72e90210697b7ced0a0d2b7559868f35db3dd Mon Sep 17 00:00:00 2001 From: mathurinm Date: Wed, 8 Dec 2021 13:32:55 +0100 Subject: [PATCH 14/34] more comments on glmnet behavior --- solvers/glmnet.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 8b058446..6fd24038 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -29,6 +29,8 @@ class Solver(BaseSolver): stopping_criterion = SufficientProgressCriterion( patience=7, eps=1e-38, strategy='tolerance') + # We use the tolerance strategy because if maxit is too low and glmnet + # convergence check fails, it returns an empty model def skip(self, X, y, lmbd, fit_intercept): # XXX - glmnet support intercept, adapt the API @@ -44,7 +46,7 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.glmnet = robjects.r['glmnet'] def run(self, tol): - # even if maxit=0, glmnet can return non zero coefficients. To get the + # Even if maxit=0, glmnet can return non zero coefficients. To get the # initial point on the curve, we set tol=0: this way, glmnet # convergence check fails, and rather than returning the current # iterate, it returns a 0 model. @@ -57,6 +59,10 @@ def run(self, tol): # plateau before convergence thresh = tol ** 1.8 + # We fit with a single lambda because when computing a path, glmnet + # early stops the path based on an uncontrollable criterion. Fitting + # with a single lambda may be suboptimal if it is small, but there is + # no other way to force glmnet to solve for a prescribed lambda. fit_dict = {"lambda": self.lmbd / len(self.y)} glmnet_fit = self.glmnet(self.X, self.y, intercept=False, From 96e529203ff800a05d6963c4d84158b0d79d00fc Mon Sep 17 00:00:00 2001 From: Thomas Moreau Date: Thu, 9 Dec 2021 10:58:04 +0100 Subject: [PATCH 15/34] Update solvers/glmnet.py --- solvers/glmnet.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 6fd24038..a11383e6 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -27,10 +27,11 @@ class Solver(BaseSolver): ] support_sparse = False - stopping_criterion = SufficientProgressCriterion( - patience=7, eps=1e-38, strategy='tolerance') # We use the tolerance strategy because if maxit is too low and glmnet # convergence check fails, it returns an empty model + stopping_criterion = SufficientProgressCriterion( + patience=7, eps=1e-38, strategy='tolerance' + ) def skip(self, X, y, lmbd, fit_intercept): # XXX - glmnet support intercept, adapt the API From ec88f6531ab3c99098357bff85708638678efade Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 11 May 2022 09:56:55 +0200 Subject: [PATCH 16/34] ENH add Lasso.jl Julia solver --- solvers/lasso_jl.jl | 36 ++++++++++++++++++++ solvers/lasso_jl.py | 80 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 solvers/lasso_jl.jl create mode 100644 solvers/lasso_jl.py diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl new file mode 100644 index 00000000..3b8e02c5 --- /dev/null +++ b/solvers/lasso_jl.jl @@ -0,0 +1,36 @@ +using Distributions +using GLM +using Lasso +using LinearAlgebra +using PyCall +using SparseArrays + +function scipyCSC_to_julia(A) + m, n = A.shape + colPtr = Int[i + 1 for i in PyArray(A."indptr")] + rowVal = Int[i + 1 for i in PyArray(A."indices")] + nzVal = Vector{Float64}(PyArray(A."data")) + B = SparseMatrixCSC{Float64,Int}(m, n, colPtr, rowVal, nzVal) + + return B +end + +function solve_lasso( + X, y::Vector{Float64}, lambda::Vector{Float64}, fit_intercept::Bool, tol::Float64 +) + lasso_fit = fit( + LassoPath, + X, + y; + λ=lambda, + standardize=false, + intercept=fit_intercept, + randomize=false, + stopearly=false, + maxncoef=max(size(X, 1), size(X, 2)) * 100, + cd_tol=tol, + ) + w = coef(lasso_fit) + + return w +end diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py new file mode 100644 index 00000000..ab3d0f5d --- /dev/null +++ b/solvers/lasso_jl.py @@ -0,0 +1,80 @@ +from pathlib import Path + +import numpy as np +from benchopt import safe_import_context +from benchopt.helpers.julia import (JuliaSolver, assert_julia_installed, + get_jl_interpreter) +from benchopt.runner import INFINITY +from benchopt.stopping_criterion import SufficientProgressCriterion + +with safe_import_context() as import_ctx: + assert_julia_installed() + from scipy import sparse + + +# File containing the function to be called from julia +JULIA_SOLVER_FILE = str(Path(__file__).with_suffix(".jl")) + + +class Solver(JuliaSolver): + name = "lasso_jl" + stopping_criterion = SufficientProgressCriterion( + patience=7, eps=1e-15, strategy="tolerance" + ) + julia_requirements = [ + "Distributions", + "GLM", + "Lasso", + "LinearAlgebra", + "PyCall", + "SparseArrays", + ] + references = [ + 'J. Friedman, T. J. Hastie and R. Tibshirani, "Regularization paths ' + 'for generalized linear models via coordinate descent", ' + "J. Stat. Softw., vol. 33, no. 1, pp. 1-22, NIH Public Access (2010)" + ] + support_sparse = True + + def set_objective(self, X, y, lmbd, fit_intercept): + self.n, self.p = X.shape + self.X = X + self.y = y + self.lmbd = np.array([lmbd]) + self.fit_intercept = fit_intercept + + jl = get_jl_interpreter() + jl.include(str(JULIA_SOLVER_FILE)) + self.solve_lasso = jl.solve_lasso + + if sparse.issparse(X): + scipyCSC_to_julia = jl.pyfunctionret( + jl.scipyCSC_to_julia, jl.Any, jl.PyObject + ) + self.X = scipyCSC_to_julia(X) + + # Trigger Julia JIT compilation + self.run(1e-5) + + def run(self, tol): + if tol == INFINITY: + # Lasso.jl always runs one iteration before checking convergence, + # so we use this workaround to force the solver to start at zero + coefs_dim = self.p + 1 if self.fit_intercept else self.p + self.coefs = np.zeros(coefs_dim) + else: + self.coefs = self.solve_lasso( + self.X, + self.y, + self.lmbd / len(self.y), + self.fit_intercept, + tol**1.8, + ) + + def get_result(self): + coefs = np.ravel(self.coefs) + + if self.fit_intercept: + coefs = np.hstack((coefs[1:], coefs[0])) + + return coefs From bf091865f0bbff52ec21e3371d54e54d053a677a Mon Sep 17 00:00:00 2001 From: Johan Larsson <13087841+jolars@users.noreply.github.com> Date: Wed, 11 May 2022 14:37:58 +0200 Subject: [PATCH 17/34] FIX remove unnecessary using calls Co-authored-by: Thomas Moreau --- solvers/lasso_jl.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl index 3b8e02c5..9f896906 100644 --- a/solvers/lasso_jl.jl +++ b/solvers/lasso_jl.jl @@ -1,7 +1,4 @@ -using Distributions -using GLM using Lasso -using LinearAlgebra using PyCall using SparseArrays From dc0d5be5993e1fa40db5bc10275469da29b6404b Mon Sep 17 00:00:00 2001 From: Johan Larsson <13087841+jolars@users.noreply.github.com> Date: Wed, 11 May 2022 14:38:09 +0200 Subject: [PATCH 18/34] FIX remove unnecessary dependencies Co-authored-by: Thomas Moreau --- solvers/lasso_jl.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index ab3d0f5d..979e21e9 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -22,10 +22,7 @@ class Solver(JuliaSolver): patience=7, eps=1e-15, strategy="tolerance" ) julia_requirements = [ - "Distributions", - "GLM", "Lasso", - "LinearAlgebra", "PyCall", "SparseArrays", ] From b60f9057c1d4dcdcb95e6ea8411adefbedae1d37 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 11 May 2022 15:22:26 +0200 Subject: [PATCH 19/34] fix: set randomize and stopearly to defaults --- solvers/lasso_jl.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl index 3b8e02c5..0b93913a 100644 --- a/solvers/lasso_jl.jl +++ b/solvers/lasso_jl.jl @@ -25,8 +25,6 @@ function solve_lasso( λ=lambda, standardize=false, intercept=fit_intercept, - randomize=false, - stopearly=false, maxncoef=max(size(X, 1), size(X, 2)) * 100, cd_tol=tol, ) From 406b0718239b8755d374ca57fda1e915ca3c628b Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 11 May 2022 15:48:55 +0200 Subject: [PATCH 20/34] fix: run tol == INFINITY iteration on julia side --- solvers/lasso_jl.jl | 35 +++++++++++++++++++++++------------ solvers/lasso_jl.py | 21 ++++++++------------- 2 files changed, 31 insertions(+), 25 deletions(-) diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl index 8113edb2..11c497b7 100644 --- a/solvers/lasso_jl.jl +++ b/solvers/lasso_jl.jl @@ -13,19 +13,30 @@ function scipyCSC_to_julia(A) end function solve_lasso( - X, y::Vector{Float64}, lambda::Vector{Float64}, fit_intercept::Bool, tol::Float64 + X, + y::Vector{Float64}, + lambda::Vector{Float64}, + fit_intercept::Bool, + tol::Float64, + get_null_solution::Bool, ) - lasso_fit = fit( - LassoPath, - X, - y; - λ=lambda, - standardize=false, - intercept=fit_intercept, - maxncoef=max(size(X, 1), size(X, 2)) * 100, - cd_tol=tol, - ) - w = coef(lasso_fit) + p = size(X, 2) + + w = if fit_intercept zeros(Float64, p + 1) else zeros(Float64, p) end + + if !get_null_solution + lasso_fit = fit( + LassoPath, + X, + y; + λ=lambda, + standardize=false, + intercept=fit_intercept, + maxncoef=max(size(X, 1), size(X, 2)) * 100, + cd_tol=tol, + ) + w = coef(lasso_fit) + end return w end diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 979e21e9..7f3d5e93 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -54,19 +54,14 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.run(1e-5) def run(self, tol): - if tol == INFINITY: - # Lasso.jl always runs one iteration before checking convergence, - # so we use this workaround to force the solver to start at zero - coefs_dim = self.p + 1 if self.fit_intercept else self.p - self.coefs = np.zeros(coefs_dim) - else: - self.coefs = self.solve_lasso( - self.X, - self.y, - self.lmbd / len(self.y), - self.fit_intercept, - tol**1.8, - ) + self.coefs = self.solve_lasso( + self.X, + self.y, + self.lmbd / len(self.y), + self.fit_intercept, + tol**1.8, + tol == INFINITY + ) def get_result(self): coefs = np.ravel(self.coefs) From b48933fded6ae5475f8d52574fb82283d4ed7289 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 11 May 2022 15:49:29 +0200 Subject: [PATCH 21/34] fix: decrease tolerance used for JIT compilation --- solvers/lasso_jl.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 7f3d5e93..1547e026 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -51,7 +51,7 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.X = scipyCSC_to_julia(X) # Trigger Julia JIT compilation - self.run(1e-5) + self.run(1e-2) def run(self, tol): self.coefs = self.solve_lasso( From 77ac79f9a56c8080dba71acf0d1e7b056b24af46 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 11:36:29 +0200 Subject: [PATCH 22/34] fix: catch convergence errors and return prev sol Lasso.jl throws an error when the algorithm fails to converge and since maximum iterations cannot be controlled, we hav to catch the error and return the previous solution instead. --- solvers/lasso_jl.jl | 31 +++++++++++++++++++------------ solvers/lasso_jl.py | 11 +++++++++-- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl index 11c497b7..0ef2b64e 100644 --- a/solvers/lasso_jl.jl +++ b/solvers/lasso_jl.jl @@ -23,20 +23,27 @@ function solve_lasso( p = size(X, 2) w = if fit_intercept zeros(Float64, p + 1) else zeros(Float64, p) end + + converged = true + # TODO(jolars): once https://github.com/JuliaStats/Lasso.jl/issues/70 or + # maybe https://github.com/JuliaStats/Lasso.jl/issues/71 is/are resolved, + # we should not need the try-catch here if !get_null_solution - lasso_fit = fit( - LassoPath, - X, - y; - λ=lambda, - standardize=false, - intercept=fit_intercept, - maxncoef=max(size(X, 1), size(X, 2)) * 100, - cd_tol=tol, - ) - w = coef(lasso_fit) + try + lasso_fit = fit( + LassoPath, + X, + y; + λ=lambda, standardize=false, intercept=fit_intercept, + maxncoef=max(size(X, 1), size(X, 2)) * 100, + cd_tol=tol, + ) + w = coef(lasso_fit) + catch + converged = false + end end - return w + return w, converged end diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 1547e026..3eaf2100 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -51,17 +51,24 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.X = scipyCSC_to_julia(X) # Trigger Julia JIT compilation + self.prev_solution = np.zeros(self.p + 1) if fit_intercept else np.zeros(self.p) self.run(1e-2) def run(self, tol): - self.coefs = self.solve_lasso( + coefs, converged = self.solve_lasso( self.X, self.y, self.lmbd / len(self.y), self.fit_intercept, tol**1.8, - tol == INFINITY + tol == INFINITY, ) + if converged: + self.coefs = coefs + self.prev_solution = coefs + else: + self.coefs = self.prev_solution + def get_result(self): coefs = np.ravel(self.coefs) From 04d400beff0aa66fe62946f63c2bb0c620acc1b5 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 12:14:49 +0200 Subject: [PATCH 23/34] fix(style): resolve linter warnings --- solvers/lasso_jl.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 3eaf2100..1598b14a 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -2,8 +2,11 @@ import numpy as np from benchopt import safe_import_context -from benchopt.helpers.julia import (JuliaSolver, assert_julia_installed, - get_jl_interpreter) +from benchopt.helpers.julia import ( + JuliaSolver, + assert_julia_installed, + get_jl_interpreter, +) from benchopt.runner import INFINITY from benchopt.stopping_criterion import SufficientProgressCriterion @@ -51,7 +54,8 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.X = scipyCSC_to_julia(X) # Trigger Julia JIT compilation - self.prev_solution = np.zeros(self.p + 1) if fit_intercept else np.zeros(self.p) + w_dim = self.p + 1 if fit_intercept else self.p + self.prev_solution = np.zeros(w_dim) self.run(1e-2) def run(self, tol): @@ -69,7 +73,6 @@ def run(self, tol): else: self.coefs = self.prev_solution - def get_result(self): coefs = np.ravel(self.coefs) From f05cfc8f74fb52ce0afd6bb1567086e2de4df82b Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 12:24:22 +0200 Subject: [PATCH 24/34] fix(ci): fix julia error by disabling multithreads --- .github/workflows/main.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b8ac4e88..5c9bfbd5 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -62,6 +62,11 @@ jobs: branch=${BENCHOPT_BRANCH##*:} pip install -U git+https://github.com/$user/benchopt@$branch + - name: Disable Julia multi-threading on MacOS + if: matrix.os == 'macos-latest' + env: + JULIA_NUM_THREADS: 1 + - name: Test run: | export OMP_NUM_THREADS=1 From 19b99cb69ba812b19f20b5a328349126eb99352c Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 12:33:20 +0200 Subject: [PATCH 25/34] Revert "fix(ci): fix julia error by disabling multithreads" This reverts commit f05cfc8f74fb52ce0afd6bb1567086e2de4df82b. --- .github/workflows/main.yml | 5 ----- 1 file changed, 5 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5c9bfbd5..b8ac4e88 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -62,11 +62,6 @@ jobs: branch=${BENCHOPT_BRANCH##*:} pip install -U git+https://github.com/$user/benchopt@$branch - - name: Disable Julia multi-threading on MacOS - if: matrix.os == 'macos-latest' - env: - JULIA_NUM_THREADS: 1 - - name: Test run: | export OMP_NUM_THREADS=1 From 704f1eae51b7943af0545d7b2a315c5364276c54 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 13:53:27 +0200 Subject: [PATCH 26/34] fix(ci): try to fix failing MacOS julia tests --- .github/workflows/main.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b8ac4e88..55720bb8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -62,6 +62,11 @@ jobs: branch=${BENCHOPT_BRANCH##*:} pip install -U git+https://github.com/$user/benchopt@$branch + - name: Disable Julia multi-threading on MacOS + if: matrix.os == 'macos-latest' + run: | + export JULIA_NUM_THREADS=1 + - name: Test run: | export OMP_NUM_THREADS=1 From b4f3c2aa624ef7de80a38ace5b4ad22c55ce22a0 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 14:00:23 +0200 Subject: [PATCH 27/34] fix(ci): disable julia multi-threading on osx --- .github/workflows/main.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b8ac4e88..55720bb8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -62,6 +62,11 @@ jobs: branch=${BENCHOPT_BRANCH##*:} pip install -U git+https://github.com/$user/benchopt@$branch + - name: Disable Julia multi-threading on MacOS + if: matrix.os == 'macos-latest' + run: | + export JULIA_NUM_THREADS=1 + - name: Test run: | export OMP_NUM_THREADS=1 From 1a2f1f6e84068f16ca50275b0076b5a6a3814bb6 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 12 May 2022 14:15:31 +0200 Subject: [PATCH 28/34] fix(ci): move JULIA_NUM_THREADS export to Test run --- .github/workflows/main.yml | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 55720bb8..75ce95a7 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -62,14 +62,10 @@ jobs: branch=${BENCHOPT_BRANCH##*:} pip install -U git+https://github.com/$user/benchopt@$branch - - name: Disable Julia multi-threading on MacOS - if: matrix.os == 'macos-latest' - run: | - export JULIA_NUM_THREADS=1 - - name: Test run: | export OMP_NUM_THREADS=1 + export JULIA_NUM_THREADS=1 benchopt test . --env-name bench_test_env -vl benchopt test . --env-name bench_test_env -vl --skip-install From a5346f101099c13206eca7138b4142f522545501 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Fri, 13 May 2022 10:07:04 +0200 Subject: [PATCH 29/34] fix(ci): stop testing on OSX --- test_config.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test_config.py b/test_config.py index baf49c4a..a2f1f5b8 100644 --- a/test_config.py +++ b/test_config.py @@ -13,6 +13,10 @@ def check_test_solver_install(solver_class): if 'julia' in solver_class.name.lower(): pytest.xfail('Julia install from conda fails currently.') + # Julia segfauls on OSX + if solver_class.name.lower() == 'lasso_jl' and sys.platform == 'darwin': + pytest.xfail('Julia solvers segfauls on OSX.') + # ModOpt install change numpy version, breaking celer install. # See CEA-COSMIC/ModOpt#144. Skipping for now if ('modopt' in solver_class.name.lower()): From 4f09c92667c84b16f04adc5f7d21dbcb57de47f7 Mon Sep 17 00:00:00 2001 From: mathurinm Date: Tue, 17 May 2022 11:11:30 +0200 Subject: [PATCH 30/34] imports in import_ctx --- solvers/lasso_jl.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 1598b14a..c8f30b49 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -1,6 +1,3 @@ -from pathlib import Path - -import numpy as np from benchopt import safe_import_context from benchopt.helpers.julia import ( JuliaSolver, @@ -12,11 +9,12 @@ with safe_import_context() as import_ctx: assert_julia_installed() + import numpy as np from scipy import sparse - -# File containing the function to be called from julia -JULIA_SOLVER_FILE = str(Path(__file__).with_suffix(".jl")) + from pathlib import Path + # File containing the function to be called from julia + JULIA_SOLVER_FILE = str(Path(__file__).with_suffix(".jl")) class Solver(JuliaSolver): @@ -34,7 +32,6 @@ class Solver(JuliaSolver): 'for generalized linear models via coordinate descent", ' "J. Stat. Softw., vol. 33, no. 1, pp. 1-22, NIH Public Access (2010)" ] - support_sparse = True def set_objective(self, X, y, lmbd, fit_intercept): self.n, self.p = X.shape From b50db118fc3c3881689a2218c29df79712bd9421 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 18 May 2022 12:28:01 +0200 Subject: [PATCH 31/34] FIX ignore FutureWarnings from pyjulia --- solvers/lasso_jl.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index c8f30b49..2f58e778 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -9,6 +9,7 @@ with safe_import_context() as import_ctx: assert_julia_installed() + import warnings import numpy as np from scipy import sparse @@ -53,9 +54,16 @@ def set_objective(self, X, y, lmbd, fit_intercept): # Trigger Julia JIT compilation w_dim = self.p + 1 if fit_intercept else self.p self.prev_solution = np.zeros(w_dim) + + warnings.filterwarnings("ignore", category=FutureWarning) self.run(1e-2) def run(self, tol): + # remove possibly spurious warnings from pyjulia + # TODO: remove filter when https://github.com/JuliaPy/pyjulia/issues/497 + # is resolved or otherwise fix the warning + warnings.filterwarnings("ignore", category=FutureWarning) + coefs, converged = self.solve_lasso( self.X, self.y, From 9661be3fd5bd82979a28477fe18083f5069af696 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 18 May 2022 12:29:59 +0200 Subject: [PATCH 32/34] FIX shorten line to appease flake8 --- solvers/lasso_jl.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 2f58e778..dfe78201 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -60,8 +60,9 @@ def set_objective(self, X, y, lmbd, fit_intercept): def run(self, tol): # remove possibly spurious warnings from pyjulia - # TODO: remove filter when https://github.com/JuliaPy/pyjulia/issues/497 - # is resolved or otherwise fix the warning + # TODO: remove filter when + # https://github.com/JuliaPy/pyjulia/issues/497 is resolved or + # otherwise fix the warning warnings.filterwarnings("ignore", category=FutureWarning) coefs, converged = self.solve_lasso( From fad0d0c31c045593740d0edbe82a738b3c9e89ff Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Wed, 18 May 2022 12:31:23 +0200 Subject: [PATCH 33/34] FIX remove trailing whitespace --- solvers/lasso_jl.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index dfe78201..453af5d4 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -61,7 +61,7 @@ def set_objective(self, X, y, lmbd, fit_intercept): def run(self, tol): # remove possibly spurious warnings from pyjulia # TODO: remove filter when - # https://github.com/JuliaPy/pyjulia/issues/497 is resolved or + # https://github.com/JuliaPy/pyjulia/issues/497 is resolved or # otherwise fix the warning warnings.filterwarnings("ignore", category=FutureWarning) From ed18ec998b5e893bd474074b7c1b908b18941c27 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Mon, 20 Jun 2022 16:44:29 +0200 Subject: [PATCH 34/34] refactor: increase the maximum number of iters --- solvers/lasso_jl.jl | 2 ++ solvers/lasso_jl.py | 1 + 2 files changed, 3 insertions(+) diff --git a/solvers/lasso_jl.jl b/solvers/lasso_jl.jl index 0ef2b64e..715f52a7 100644 --- a/solvers/lasso_jl.jl +++ b/solvers/lasso_jl.jl @@ -18,6 +18,7 @@ function solve_lasso( lambda::Vector{Float64}, fit_intercept::Bool, tol::Float64, + cd_maxiter::Int, get_null_solution::Bool, ) p = size(X, 2) @@ -37,6 +38,7 @@ function solve_lasso( y; λ=lambda, standardize=false, intercept=fit_intercept, maxncoef=max(size(X, 1), size(X, 2)) * 100, + cd_maxiter=cd_maxiter, cd_tol=tol, ) w = coef(lasso_fit) diff --git a/solvers/lasso_jl.py b/solvers/lasso_jl.py index 453af5d4..218d176d 100644 --- a/solvers/lasso_jl.py +++ b/solvers/lasso_jl.py @@ -71,6 +71,7 @@ def run(self, tol): self.lmbd / len(self.y), self.fit_intercept, tol**1.8, + 1_000_000, tol == INFINITY, ) if converged: