diff --git a/datasets/dct_inpainting.py b/datasets/dct_inpainting.py new file mode 100644 index 0000000..fab4a27 --- /dev/null +++ b/datasets/dct_inpainting.py @@ -0,0 +1,61 @@ +from benchopt import BaseDataset +from benchopt import safe_import_context + +with safe_import_context() as import_ctx: + import pooch # noqa: F401 + import numpy as np + from scipy.datasets import ascent + from scipy.fftpack import dctn, idctn + from scipy.sparse.linalg import LinearOperator + + +class Dataset(BaseDataset): + name = "dct_inpainting" + requirements = ["pooch"] + + def __init__(self, random_state=27): + self.random_state = random_state + self.X, self.y = None, None + + @staticmethod + def _load_image(noise_level=0.01, random_state=27): + rng = np.random.RandomState(random_state) + + y = ascent().astype(np.float64) / 255.0 + y += noise_level * rng.normal(0, 1.0, size=y.shape) + + return y + + def get_data(self): + if self.X is None or self.y is None: + y = self._load_image(noise_level=0.1, + random_state=self.random_state) + pw, ph = y.shape + + rng = np.random.RandomState(self.random_state) + mask = rng.binomial(1, 0.9, pw * ph).astype(bool) + + y = y.flatten() + y[mask] = 0.0 + + def _fwd_operator(u): + dct_coeffs = idctn(u.reshape((pw, ph)), norm="ortho").reshape( + (pw * ph, )) + dct_coeffs[mask] = 0.0 + return dct_coeffs + + def _rfwd_operator(x): + masked_coeffs = x.copy() + masked_coeffs[mask] = 0.0 + idct_coeffs = dctn(masked_coeffs.reshape((pw, ph)), + norm="ortho").reshape((pw * ph, )) + return idct_coeffs + + _fwd_scipy_operator = LinearOperator( + (pw * ph, pw * ph), + matvec=_fwd_operator, + rmatvec=_rfwd_operator, + ) + + self.X, self.y = _fwd_scipy_operator, y + return dict(X=self.X, y=self.y) diff --git a/solvers/blitz.py b/solvers/blitz.py index 753900e..11ab6a8 100644 --- a/solvers/blitz.py +++ b/solvers/blitz.py @@ -3,6 +3,7 @@ with safe_import_context() as import_ctx: + from scipy.sparse.linalg import LinearOperator import numpy as np import blitzl1 @@ -19,6 +20,21 @@ class Solver(BaseSolver): 'pip:git+https://github.com/tbjohns/blitzl1.git@master' ] + references = [ + 'T. B. Johnson and C. Guestrin, "Blitz: A Principled Meta-Algorithm ' + 'for Scaling Sparse Optimization", ICML, ' + 'vol. 37, pp. 1171-1179 (2015)' + ] + + def skip(self, X, y, lmbd, fit_intercept): + if fit_intercept: + return True, f"{self.name} does not handle fit_intercept" + + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/cd.py b/solvers/cd.py index 02007cc..56208b3 100644 --- a/solvers/cd.py +++ b/solvers/cd.py @@ -3,6 +3,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from scipy import sparse from numba import njit @@ -37,6 +38,8 @@ def skip(self, X, y, lmbd, fit_intercept): # XXX - not implemented but this should be quite easy if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" return False, None diff --git a/solvers/celer.py b/solvers/celer.py index b2527ee..c9023db 100644 --- a/solvers/celer.py +++ b/solvers/celer.py @@ -6,6 +6,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from celer import Lasso from sklearn.exceptions import ConvergenceWarning @@ -20,6 +21,12 @@ class Solver(BaseSolver): install_cmd = 'conda' requirements = ['pip:celer'] + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/cuml.py b/solvers/cuml.py index a85c6da..ab5ef2f 100644 --- a/solvers/cuml.py +++ b/solvers/cuml.py @@ -5,6 +5,7 @@ cuda_version = None with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from scipy import sparse cuda_version = requires_gpu() @@ -42,6 +43,12 @@ class Solver(BaseSolver): eps=1e-12, patience=5, strategy='iteration' ) + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd if sparse.issparse(X): diff --git a/solvers/cyanure.py b/solvers/cyanure.py index 675be5b..c15b7ee 100644 --- a/solvers/cyanure.py +++ b/solvers/cyanure.py @@ -4,6 +4,7 @@ with safe_import_context() as import_ctx: import scipy + from scipy.sparse.linalg import LinearOperator import numpy as np from cyanure import Regression @@ -19,6 +20,12 @@ class Solver(BaseSolver): 'Arxiv eprint 1912.08165 (2019)' ] + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/glmnet.py b/solvers/glmnet.py index 7c9ac65..fb3ba3f 100644 --- a/solvers/glmnet.py +++ b/solvers/glmnet.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: import numpy as np from scipy import sparse + from scipy.sparse.linalg import LinearOperator from rpy2 import robjects from rpy2.robjects import numpy2ri, packages @@ -34,6 +35,12 @@ class Solver(BaseSolver): patience=7, eps=1e-38, strategy='tolerance' ) + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): if sparse.issparse(X): r_Matrix = packages.importr("Matrix") diff --git a/solvers/glum.py b/solvers/glum.py index 9c6d79d..8275ec3 100644 --- a/solvers/glum.py +++ b/solvers/glum.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from glum import GeneralizedLinearRegressor from sklearn.exceptions import ConvergenceWarning @@ -39,6 +40,9 @@ def set_objective(self, X, y, lmbd, fit_intercept): def skip(self, X, y, lmbd, fit_intercept): # glum instantiates a Hessian of size (n_features, n_features) + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + if X.shape[1] > 20_000: return True, "glum does not support n_features >= 20000" return False, None diff --git a/solvers/irls.py b/solvers/irls.py index f7b5294..2ec231b 100644 --- a/solvers/irls.py +++ b/solvers/irls.py @@ -3,6 +3,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from sklearn.linear_model import Ridge @@ -28,6 +29,9 @@ def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + return False, None def set_objective(self, X, y, lmbd, fit_intercept): diff --git a/solvers/julia_pgd.py b/solvers/julia_pgd.py index 999a419..19ccead 100644 --- a/solvers/julia_pgd.py +++ b/solvers/julia_pgd.py @@ -6,6 +6,7 @@ from benchopt.helpers.julia import assert_julia_installed with safe_import_context() as import_ctx: + from scipy.sparse.linalg import LinearOperator assert_julia_installed() @@ -33,6 +34,9 @@ def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + return False, None def set_objective(self, X, y, lmbd, fit_intercept): diff --git a/solvers/l_bfgs_b.py b/solvers/l_bfgs_b.py index 2410c35..c1b2ff5 100644 --- a/solvers/l_bfgs_b.py +++ b/solvers/l_bfgs_b.py @@ -3,6 +3,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from numpy.linalg import norm from scipy.optimize import fmin_l_bfgs_b @@ -31,6 +32,8 @@ def skip(self, X, y, lmbd, fit_intercept): # XXX - not implemented but this should be quite easy if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" return False, None diff --git a/solvers/lars.py b/solvers/lars.py index e2ad6fc..4afb339 100644 --- a/solvers/lars.py +++ b/solvers/lars.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: import warnings import numpy as np + from scipy.sparse.linalg import LinearOperator from sklearn.linear_model import LassoLars from sklearn.exceptions import ConvergenceWarning @@ -21,6 +22,12 @@ class Solver(BaseSolver): ] support_sparse = False + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/lightning.py b/solvers/lightning.py index bf7be8e..e94b979 100644 --- a/solvers/lightning.py +++ b/solvers/lightning.py @@ -4,6 +4,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from lightning.regression import CDRegressor @@ -27,6 +28,8 @@ class Solver(BaseSolver): def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" return False, None diff --git a/solvers/modopt_fista.py b/solvers/modopt_fista.py index 9443f65..cb31590 100644 --- a/solvers/modopt_fista.py +++ b/solvers/modopt_fista.py @@ -6,6 +6,8 @@ with safe_import_context() as import_ctx: from scipy import sparse + from scipy.sparse.linalg import LinearOperator + from scipy.linalg.interpolative import estimate_spectral_norm from modopt.opt.algorithms import ForwardBackward from modopt.opt.proximity import SparseThreshold from modopt.opt.linear import Identity @@ -45,7 +47,9 @@ def set_objective(self, X, y, lmbd, fit_intercept): self.var_init = np.zeros(n_features) def run(self, callback): - if sparse.issparse(self.X): + if isinstance(self.X, LinearOperator): + L = estimate_spectral_norm(self.X) ** 2 + elif sparse.issparse(self.X): L = sparse.linalg.svds(self.X, k=1)[1][0] ** 2 else: L = np.linalg.norm(self.X, ord=2) ** 2 @@ -64,6 +68,7 @@ def run(self, callback): q_lazy = 1 / 10 n_features = self.X.shape[1] + if self.fit_intercept: def op(w): return self.X @ w[:n_features] + w[-1] diff --git a/solvers/modopt_pogm.py b/solvers/modopt_pogm.py index 168db17..25bc046 100644 --- a/solvers/modopt_pogm.py +++ b/solvers/modopt_pogm.py @@ -6,6 +6,8 @@ with safe_import_context() as import_ctx: from scipy import sparse + from scipy.sparse.linalg import LinearOperator + from scipy.linalg.interpolative import estimate_spectral_norm from modopt.opt.algorithms import POGM from modopt.opt.proximity import SparseThreshold from modopt.opt.linear import Identity @@ -63,7 +65,9 @@ def trans_op(res): weights = np.full(self.X.shape[1], self.lmbd) - if sparse.issparse(self.X): + if isinstance(self.X, LinearOperator): + L = estimate_spectral_norm(self.X) ** 2 + elif sparse.issparse(self.X): L = sparse.linalg.svds(self.X, k=1)[1][0] ** 2 else: L = np.linalg.norm(self.X, ord=2) ** 2 diff --git a/solvers/noncvx_pro.py b/solvers/noncvx_pro.py index 8317a48..3bb1c5c 100644 --- a/solvers/noncvx_pro.py +++ b/solvers/noncvx_pro.py @@ -7,6 +7,7 @@ from numpy.linalg import norm import scipy.optimize as sciop from scipy.sparse import issparse + from scipy.sparse.linalg import LinearOperator class Solver(BaseSolver): @@ -30,6 +31,9 @@ def set_objective(self, X, y, lmbd, fit_intercept): def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + return False, None def run(self, n_iter): diff --git a/solvers/python_pgd.py b/solvers/python_pgd.py index 81f600b..963fd24 100644 --- a/solvers/python_pgd.py +++ b/solvers/python_pgd.py @@ -3,6 +3,8 @@ with safe_import_context() as import_ctx: import numpy as np from scipy import sparse + from scipy.linalg.interpolative import estimate_spectral_norm + from scipy.sparse.linalg import LinearOperator class Solver(BaseSolver): @@ -63,7 +65,9 @@ def get_result(self): return self.w def compute_lipschitz_constant(self): - if not sparse.issparse(self.X): + if isinstance(self.X, LinearOperator): + L = estimate_spectral_norm(self.X) ** 2 + elif not sparse.issparse(self.X): L = np.linalg.norm(self.X, ord=2) ** 2 else: L = sparse.linalg.svds(self.X, k=1)[1][0] ** 2 diff --git a/solvers/r_pgd.py b/solvers/r_pgd.py index 54c6542..531f91d 100644 --- a/solvers/r_pgd.py +++ b/solvers/r_pgd.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from rpy2 import robjects from rpy2.robjects import numpy2ri @@ -36,6 +37,8 @@ class Solver(BaseSolver): def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" return False, None diff --git a/solvers/skglm.py b/solvers/skglm.py index 53d09ce..bb0b119 100644 --- a/solvers/skglm.py +++ b/solvers/skglm.py @@ -4,6 +4,7 @@ with safe_import_context() as import_ctx: import warnings import numpy as np + from scipy.sparse.linalg import LinearOperator from skglm import Lasso from sklearn.exceptions import ConvergenceWarning @@ -21,6 +22,21 @@ class Solver(BaseSolver): 'pip:skglm' ] + references = [ + 'Q. Bertrand and Q. Klopfenstein and P.-A. Bannier and G. Gidel' + 'and M. Massias' + '"Beyond L1: Faster and Better Sparse Models with skglm", ' + 'https://arxiv.org/abs/2204.07826' + ] + + def skip(self, X, y, lmbd, fit_intercept): + if fit_intercept: + return True, f"{self.name} does not handle fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/sklearn.py b/solvers/sklearn.py index 177f4eb..a552230 100644 --- a/solvers/sklearn.py +++ b/solvers/sklearn.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: import numpy as np + from scipy.sparse.linalg import LinearOperator from sklearn.linear_model import Lasso from sklearn.exceptions import ConvergenceWarning @@ -20,6 +21,12 @@ class Solver(BaseSolver): install_cmd = 'conda' requirements = ['scikit-learn'] + def skip(self, X, y, lmbd, fit_intercept): + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + + return False, None + def set_objective(self, X, y, lmbd, fit_intercept): self.X, self.y, self.lmbd = X, y, lmbd self.fit_intercept = fit_intercept diff --git a/solvers/snapml.py b/solvers/snapml.py index 14866d2..15a3edf 100644 --- a/solvers/snapml.py +++ b/solvers/snapml.py @@ -5,6 +5,7 @@ with safe_import_context() as import_ctx: from snapml import LinearRegression import numpy as np + from scipy.sparse.linalg import LinearOperator class Solver(BaseSolver): @@ -26,6 +27,9 @@ class Solver(BaseSolver): def skip(self, X, y, lmbd, fit_intercept): if self.gpu and get_cuda_version() is None: return True, "snapml[gpu=True] needs a GPU to run" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + return False, None def set_objective(self, X, y, lmbd, fit_intercept): diff --git a/solvers/spams.py b/solvers/spams.py index 4fd080a..1a2639b 100644 --- a/solvers/spams.py +++ b/solvers/spams.py @@ -4,6 +4,7 @@ with safe_import_context() as import_ctx: import scipy import numpy as np + from scipy.sparse.linalg import LinearOperator from spams import lasso, fistaFlat @@ -19,6 +20,9 @@ class Solver(BaseSolver): def skip(self, X, y, lmbd, fit_intercept): if fit_intercept: return True, f"{self.name} does not work with fit_intercept" + if isinstance(X, LinearOperator): + return True, f"{self.name} does not handle implicit operator" + return False, None def set_objective(self, X, y, lmbd, fit_intercept):