feat: remove dependence on nlopt and use L-BFGS-B in scipy as default…

… numeric solver (abess-team#95)

* feat: use L-BFGS-B in scipy as default numeric solver

remove dependence on nlopt

* chore: debug jax[cpu]

* format code style

docs/source/autoapi/numeric_solver.rst

@@ -17,7 +17,7 @@ Functions
 
 .. autoapisummary::
 
-    skscope.numeric_solver.convex_solver_nlopt
+    skscope.numeric_solver.convex_solver_LBFGS
 
 .. autoapimodule:: skscope.numeric_solver
     :members: 

docs/source/feature/Variants.rst

@@ -168,16 +168,16 @@ where
 
 - :math:`f(\theta)` is the objective function; 
 
-All solvers in :ref:`skscope <skscope_package>` use `nlopt <https://nlopt.readthedocs.io/en/latest/>`_ as the default numeric optimization solver for this problem. 
+All solvers in :ref:`skscope <skscope_package>` use `L-BFGS-B` algorithm in `scipy.optimize <https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html>`_ as the default numeric optimization solver for this problem. 
 
 In some cases, there may be additional constraints on the intrinsic structure of :math:`\theta`, which can be formulated as a set :math:`\mathcal{C}`:
 
 .. math::
     \arg\min_{\theta \in R^s, \theta \in \mathcal{C}} f(\theta).
 
-A typical example is the Gaussian graphical model for continuous random variables, which constrains :math:`\theta` on symmetric positive-definite spaces (see this example `<../userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.html>`__). Although ``nlopt`` cannot solve this problem, ``skscope`` provides a flexible interface that allows for its replacement. Specifically, users can change the default numerical optimization solver by properly setting the ``numeric_solver`` in the solver. 
+A typical example is the Gaussian graphical model for continuous random variables, which constrains :math:`\theta` on symmetric positive-definite spaces (see this example `<../userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.html>`__). Although the default numeric solver cannot solve this problem, ``skscope`` provides a flexible interface that allows for its replacement. Specifically, users can change the default numerical optimization solver by properly setting the ``numeric_solver`` in the solver. 
 
-    > Notice that, the accepted input of ``numeric_solver`` should have the same interface as ``skscope.numeric_solver.convex_solver_nlopt``.
+    > Notice that, the accepted input of ``numeric_solver`` should have the same interface as ``skscope.numeric_solver.convex_solver_LBFGS``.
 
 
 .. code-block:: python

docs/source/gallery/GeneralizedLinearModels/Inverse-gaussian-regression.ipynb

@@ -269,6 +269,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from skscope.numeric_solver import convex_solver_LBFGS\n",
     "def convex_solver_inverse_gaussian(\n",
     "    objective_func,\n",
     "    value_and_grad,\n",
@@ -282,7 +283,7 @@
     "    m = np.min(X @ params)\n",
     "    if m <= 0.0:\n",
     "        params[0] -= m\n",
-    "    return convex_solver_nlopt(objective_func, value_and_grad, params, optim_variable_set, data)"
+    "    return convex_solver_LBFGS(objective_func, value_and_grad, params, optim_variable_set, data)"
    ]
   },
   {

docs/source/gallery/GeneralizedLinearModels/gamma-regression.ipynb

@@ -80,7 +80,7 @@
     "import numpy as np\n",
     "import jax.numpy as jnp\n",
     "from skscope import ScopeSolver, HTPSolver, GraspSolver\n",
-    "from skscope.numeric_solver import convex_solver_nlopt"
+    "from skscope.numeric_solver import convex_solver_LBFGS"
    ]
   },
   {
@@ -284,7 +284,7 @@
     "    m = np.min(X @ params)\n",
     "    if m <= 0.0:\n",
     "        params[0] -= m\n",
-    "    return convex_solver_nlopt(objective_func, value_and_grad, params, optim_variable_set, data)\n"
+    "    return convex_solver_LBFGS(objective_func, value_and_grad, params, optim_variable_set, data)\n"
    ]
   },
   {

docs/source/userguide/install.rst

@@ -72,6 +72,5 @@ Dependencies
 The current minimum dependencies to run ``skscope`` are:
 
 - ``numpy``
-- ``nlopt``
 - ``scikit-learn>=1.2.2``
-- ``"jax[cpu]"``
+- ``jax[cpu]``

pytest/test_skmodel.py

@@ -45,7 +45,7 @@ def test_PortfolioSelection():
     # fit and test
     port = port.fit(X_train)
     score = port.score(X_test)
-    assert score > 0.05
+    assert score > 0.049
 
     # gridsearch with time-series splitting
     tscv = TimeSeriesSplit(n_splits=5)

setup.py

@@ -100,9 +100,7 @@ def build_extension(self, ext):
 
             # Multi-config generators have a different way to specify configs
             if not single_config:
-                cmake_args += [
-                    f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}"
-                ]
+                cmake_args += [f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}"]
                 build_args += ["--config", cfg]
 
         if sys.platform.startswith("darwin"):
@@ -151,7 +149,6 @@ def build_extension(self, ext):
         "numpy",
         "scikit-learn>=1.2.2",
         "jax[cpu]",
-        "nlopt",
         "scipy",
     ],
     license="MIT",

skscope/__init__.py

@@ -18,7 +18,7 @@
     OMPSolver,
 )
 from .base_solver import BaseSolver
-from .numeric_solver import convex_solver_nlopt
+from .numeric_solver import convex_solver_LBFGS
 from .skmodel import PortfolioSelection, NonlinearSelection, RobustRegression
 
 __all__ = [
@@ -30,7 +30,7 @@
     "FobaSolver",
     "ForwardSolver",
     "OMPSolver",
-    "convex_solver_nlopt",
+    "convex_solver_LBFGS",
     "PortfolioSelection",
     "NonlinearSelection",
     "RobustRegression",

skscope/base_solver.py

@@ -8,7 +8,7 @@
 from sklearn.model_selection import KFold
 import numpy as np
 import jax
-from .numeric_solver import convex_solver_nlopt
+from .numeric_solver import convex_solver_LBFGS
 import math
 from . import utilities
 
@@ -32,7 +32,7 @@ class BaseSolver(BaseEstimator):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``BaseSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -78,7 +78,7 @@ def __init__(
         sample_size=1,
         *,
         preselect=[],
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=100,
         group=None,
         ic_method=None,

skscope/layer.py

@@ -93,7 +93,7 @@ def __init__(self, dimensionality, coef=None):
     @jax.jit
     def transform_params(self, params):
         x = jnp.dot(params, self.coef)
-        return params / jnp.where(x == 0.0, 1.0, x)
+        return params / jnp.where(jnp.abs(x) < 1e-5, 1.0, x)
 
     def tree_flatten(self):
         children = (self.coef,)
@@ -133,7 +133,7 @@ def __init__(self, dimensionality, coef=None):
     def transform_params(self, params):
         p = jnp.abs(params)
         x = jnp.dot(p, self.coef)
-        return p / jnp.where(x == 0.0, 1.0, x)
+        return p / jnp.where(jnp.abs(x) < 1e-5, 1.0, x)
 
     def tree_flatten(self):
         children = (self.coef,)

skscope/numeric_solver.py

@@ -6,73 +6,31 @@
 
 import numpy as np
 import math
-import nlopt
 from scipy.optimize import minimize
 
 
-def convex_solver_nlopt(
+def convex_solver_BFGS(
     objective_func,
     value_and_grad,
     init_params,
     optim_variable_set,
     data,
 ):
-    """
-    A wrapper of ``nlopt`` solver for convex optimization.
-
-    Parameters
-    ----------
-    objective_func: callable
-        The objective function.
-        ``objective_func(params, data) -> loss``, where ``params`` is a 1-D array with shape (dimensionality,).
-    value_and_grad: callable
-        The function to compute the loss and gradient.
-        ``value_and_grad(params, data) -> (loss, grad)``, where ``params`` is a 1-D array with shape (dimensionality,).
-    init_params: array of shape (dimensionality,)
-        The initial value of the parameters to be optimized.
-    optim_variable_set: array of int
-        The index of variables to be optimized, others are fixed to the initial value.
-    data:
-        The data passed to objective_func and value_and_grad.
-
-    Returns
-    -------
-    loss: float
-        The loss of the optimized parameters, i.e., `objective_func(params, data)`.
-    optimized_params: array of shape (dimensionality,)
-        The optimized parameters.
-    """
-    best_loss = math.inf
-    best_params = None
-
-    def cache_opt_fn(x, grad):
-        nonlocal best_loss, best_params
-        # update the nonlocal variable: params
+    def fun(x):
         init_params[optim_variable_set] = x
-        if grad.size > 0:
-            loss, full_grad = value_and_grad(init_params, data)
-            grad[:] = full_grad[optim_variable_set]
-        else:
-            loss = objective_func(init_params, data)
-        if loss < best_loss:
-            best_loss = loss
-            best_params = np.copy(x)
-        return loss
+        return objective_func(init_params, data)
 
-    nlopt_solver = nlopt.opt(nlopt.LD_LBFGS, optim_variable_set.size)
-    nlopt_solver.set_min_objective(cache_opt_fn)
+    def jac(x):
+        init_params[optim_variable_set] = x
+        _, grad = value_and_grad(init_params, data)
+        return grad[optim_variable_set]
 
-    try:
-        init_params[optim_variable_set] = nlopt_solver.optimize(
-            init_params[optim_variable_set]
-        )
-        return nlopt_solver.last_optimum_value(), init_params
-    except RuntimeError:
-        init_params[optim_variable_set] = best_params
-        return best_loss, init_params
+    res = minimize(fun, init_params[optim_variable_set], method="BFGS", jac=jac)
+    init_params[optim_variable_set] = res.x
+    return res.fun, init_params
 
 
-def convex_solver_BFGS(
+def convex_solver_LBFGS(
     objective_func,
     value_and_grad,
     init_params,
@@ -88,6 +46,6 @@ def jac(x):
         _, grad = value_and_grad(init_params, data)
         return grad[optim_variable_set]
 
-    res = minimize(fun, init_params[optim_variable_set], method="BFGS", jac=jac)
+    res = minimize(fun, init_params[optim_variable_set], method="L-BFGS-B", jac=jac)
     init_params[optim_variable_set] = res.x
     return res.fun, init_params

skscope/solver.py

@@ -11,7 +11,7 @@
 import jax
 from jax import numpy as jnp
 from . import _scope, utilities
-from .numeric_solver import convex_solver_nlopt
+from .numeric_solver import convex_solver_LBFGS
 
 
 class ScopeSolver(BaseEstimator):
@@ -33,7 +33,7 @@ class ScopeSolver(BaseEstimator):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``ScopeSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=20
         Maximum number of iterations taken for converging.
     ic_method : callable, optional
@@ -122,7 +122,7 @@ def __init__(
         sample_size=1,
         *,
         preselect=[],
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=20,
         ic_method=None,
         cv=1,
@@ -699,7 +699,7 @@ class HTPSolver(BaseSolver):
         Step size of gradient descent.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``HTPSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -750,7 +750,7 @@ def __init__(
         *,
         preselect=[],
         step_size=0.005,
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=100,
         group=None,
         ic_method=None,
@@ -860,7 +860,7 @@ class IHTSolver(HTPSolver):
         Step size of gradient descent.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``IHTSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -980,7 +980,7 @@ class GraspSolver(BaseSolver):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``GraspSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -1141,7 +1141,7 @@ class FobaSolver(BaseSolver):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``FobaSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -1194,7 +1194,7 @@ def __init__(
         foba_threshold_ratio=0.5,
         strict_sparsity=True,
         preselect=[],
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=100,
         group=None,
         ic_method=None,
@@ -1416,7 +1416,7 @@ class ForwardSolver(FobaSolver):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``ForwardSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -1464,7 +1464,7 @@ def __init__(
         threshold=0.0,
         strict_sparsity=True,
         preselect=[],
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=100,
         group=None,
         ic_method=None,
@@ -1564,7 +1564,7 @@ class OMPSolver(ForwardSolver):
         An array contains the indexes of variables which must be selected.
     numeric_solver : callable, optional
         A solver for the convex optimization problem. ``OMPSolver`` will call this function to solve the convex optimization problem in each iteration.
-        It should have the same interface as ``skscope.convex_solver_nlopt``.
+        It should have the same interface as ``skscope.convex_solver_LBFGS``.
     max_iter : int, default=100
         Maximum number of iterations taken for converging.
     group : array of shape (dimensionality,), default=range(dimensionality)
@@ -1615,7 +1615,7 @@ def __init__(
         threshold=0.0,
         strict_sparsity=True,
         preselect=[],
-        numeric_solver=convex_solver_nlopt,
+        numeric_solver=convex_solver_LBFGS,
         max_iter=100,
         group=None,
         ic_method=None,