Image DCT inpainting dataset

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)

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

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
 

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

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):

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

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")

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

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):

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):

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
 

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

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
 

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]

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

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):

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