Add ARD support to bayes_search and related policies

docs/sphinx/manual/en/source/notebook/index.rst

@@ -9,6 +9,7 @@ Here, the usage of PHYSBO is introduced through tutorials.
 
    tutorial_basic
    tutorial_Gaussian_process
+   tutorial_ard
    tutorial_interactive_mode
    tutorial_once_mode
    tutorial_multi_probe

docs/sphinx/manual/ja/source/notebook/index.rst

@@ -9,6 +9,7 @@
 
    tutorial_basic
    tutorial_Gaussian_process
+   tutorial_ard
    tutorial_interactive_mode
    tutorial_once_mode
    tutorial_multi_probe

examples/single_objective/ard_example.py

@@ -0,0 +1,118 @@
+# SPDX-License-Identifier: MPL-2.0
+# Copyright (C) 2020- The University of Tokyo
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+"""
+Example: When ARD (Automatic Relevance Determination) is useful
+
+ARD is effective when the objective function depends on only a subset of
+input dimensions. The GP kernel learns a separate length scale per dimension;
+relevant dimensions get smaller length scales, irrelevant ones get larger
+(and are effectively down-weighted).
+"""
+
+import numpy as np
+import physbo
+
+weights = np.array([5.0, 1.0, 0.0, 0.0])
+D = len(weights)
+N = 1000
+M = 0
+np.random.seed(137)
+test_X = np.random.randn(N, D)
+test_X[0, :] = 0.0
+
+
+def simulator(actions: np.ndarray) -> np.ndarray:
+    """Objective function: f = -Σ_i w_i x_i^2."""
+    X2 = test_X[actions, :] ** 2
+    return -np.einsum("ai,i -> a", X2, weights)
+
+
+# ---------- Run with ARD=True ----------
+n_initial = 20
+n_bayes = 30
+n_perm = 20
+score = "EI"
+seed = 31415
+
+print("=" * 60)
+print("Running with ard=True")
+print("=" * 60)
+policy_ard = physbo.search.discrete.Policy(test_X)
+policy_ard.set_seed(seed)
+policy_ard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False)
+policy_ard.bayes_search(
+    max_num_probes=n_bayes,
+    simulator=simulator,
+    score=score,
+    ard=True,
+    is_disp=False,
+    num_rand_basis=M,
+)
+
+best_fx_ard, best_actions_ard = policy_ard.history.export_sequence_best_fx()
+best_x_ard = test_X[best_actions_ard[-1], :]
+print("\nBest value (ard=True):", best_fx_ard[-1])
+print("Best point:", best_x_ard)
+
+# ---------- Run with ARD=False ----------
+print("\n" + "=" * 60)
+print("Running with ard=False")
+print("=" * 60)
+policy_noard = physbo.search.discrete.Policy(test_X)
+policy_noard.set_seed(seed)
+policy_noard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False)
+policy_noard.bayes_search(
+    max_num_probes=n_bayes,
+    simulator=simulator,
+    score=score,
+    ard=False,
+    is_disp=False,
+)
+
+best_fx_noard, best_actions_noard = policy_noard.history.export_sequence_best_fx()
+best_x_noard = test_X[best_actions_noard[-1], :]
+print("\nBest value (ard=False):", best_fx_noard[-1])
+print("Best point:", best_x_noard)
+
+# ---------- Compare best values ----------
+print("\n" + "=" * 60)
+print("Comparison: ard=True vs ard=False")
+print("=" * 60)
+print(f"  Best value (ard=True):  {best_fx_ard[-1]:.6f}")
+print(f"  Best value (ard=False): {best_fx_noard[-1]:.6f}")
+print(f"  (Higher is better; weights are {weights}).")
+
+# ---------- Compare kernel length scales ----------
+print("\n--- Kernel length scale ---")
+ls_ard = policy_ard.get_kernel_length_scale()
+ls_noard = policy_noard.get_kernel_length_scale()
+if ls_ard is not None:
+    print("  ard=True:  per-dimension length scale (smaller = more relevant)")
+    for i in range(len(ls_ard)):
+        rel = "relevant" if i < 2 else "irrelevant"
+        print(f"    dim {i} ({rel}): {ls_ard[i]:.4f}")
+else:
+    print("  ard=True:  (not available)")
+if ls_noard is not None:
+    print("  ard=False: single length scale (isotropic kernel):", ls_noard[0])
+else:
+    print("  ard=False: (not available)")
+
+# ---------- Permutation importance (both) ----------
+print("\n--- Permutation importance ---")
+imp_mean_ard, imp_std_ard = policy_ard.get_permutation_importance(n_perm=n_perm)
+imp_mean_noard, imp_std_noard = policy_noard.get_permutation_importance(n_perm=n_perm)
+print("  ard=True:")
+for i in range(D):
+    rel = "relevant" if i < 2 else "irrelevant"
+    print(f"    dim {i} ({rel}): mean = {imp_mean_ard[i]:.4f}, std = {imp_std_ard[i]:.4f}")
+print("  ard=False:")
+for i in range(D):
+    rel = "relevant" if i < 2 else "irrelevant"
+    print(f"    dim {i} ({rel}): mean = {imp_mean_noard[i]:.4f}, std = {imp_std_noard[i]:.4f}")
+print("  (Higher mean => more important; both use the same GP predictor after optimization.)")

src/physbo/gp/core/_model.py

@@ -8,7 +8,7 @@
 import numpy as np
 
 from ... import blm
-from .. import inf
+from .. import cov, inf, lik, mean
 from . import learning
 from ._prior import Prior
 
@@ -33,6 +33,35 @@ def __init__(self, lik, mean, cov, inf="exact"):
         self.params = self.cat_params(self.lik.params, self.prior.params)
         self.stats = ()
 
+    @classmethod
+    def create_default(cls, ard=False, num_dim=None):
+        """
+        Create a default GP model with Gaussian kernel, constant mean, and Gaussian likelihood.
+
+        Parameters
+        ----------
+        ard : bool
+            If True, use ARD (per-dimension length scale) for the Gaussian kernel.
+        num_dim : int or None
+            Input dimension. Required when ard is True; optional otherwise.
+
+        Returns
+        -------
+        Model
+            A new Model instance.
+        """
+        if ard and num_dim is None:
+            raise ValueError(
+                "ard=True requires num_dim. "
+                "Provide the input dimension (e.g. from training data)."
+            )
+        cov_kernel = cov.Gauss(num_dim=num_dim, ard=ard)
+        return cls(
+            lik=lik.Gauss(),
+            mean=mean.Const(),
+            cov=cov_kernel,
+        )
+
     def cat_params(self, lik_params, prior_params):
         """
         Concatinate the likelihood and prior parameters

src/physbo/search/discrete/_policy.py

@@ -15,6 +15,7 @@
 from ._history import History
 from .. import utility
 from .. import score as search_score
+from ... import gp
 from ...gp import Predictor as gp_predictor
 from ...blm import Predictor as blm_predictor
 from ...misc import SetConfig
@@ -105,6 +106,7 @@ def __init__(self, test_X, config=None, initial_data=None, comm=None):
         else:
             self.config = config
 
+        self.ard = False
         if initial_data is not None:
             if len(initial_data) != 2:
                 msg = "ERROR: initial_data should be 2-elements tuple or list (actions and objectives)"
@@ -311,6 +313,7 @@ def bayes_search(
         num_rand_basis=0,
         optimizer=None,
         unify_method=None,
+        ard=False,
     ):
         """
         Performing Bayesian optimization.
@@ -327,7 +330,7 @@ def bayes_search(
             Base class is defined in physbo.predictor.
             If None, blm_predictor is defined.
         is_disp: bool
-             If true, process messages are outputted.
+            If true, process messages are outputted.
         simulator: callable
             Callable (function or object with ``__call__``)
             Here, action is an integer which represents the index of the candidate.
@@ -343,6 +346,9 @@ def bayes_search(
         optimizer: optimizer object, optional
             This is for compatibility with the range-based Policies.
             This is not used.
+        ard: bool
+            If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel.
+            Default is False.
 
         Returns
         -------
@@ -360,6 +366,7 @@ def bayes_search(
             simulator = None
 
         is_rand_expans = num_rand_basis != 0
+        self.ard = ard
 
         if training is not None:
             self.training = training
@@ -450,7 +457,7 @@ def get_post_fmean(self, xs):
         X = self._make_variable_X(xs)
         if self.predictor is None:
             self._warn_no_predictor("get_post_fmean()")
-            predictor = gp_predictor(self.config)
+            predictor = self._make_gp_predictor()
             predictor.fit(self.training, 0, comm=self.mpicomm)
             predictor.prepare(self.training)
             return predictor.get_post_fmean(self.training, X)
@@ -479,14 +486,55 @@ def get_post_fcov(self, xs, diag=True):
         X = self._make_variable_X(xs)
         if self.predictor is None:
             self._warn_no_predictor("get_post_fcov()")
-            predictor = gp_predictor(self.config)
+            predictor = self._make_gp_predictor()
             predictor.fit(self.training, 0, comm=self.mpicomm)
             predictor.prepare(self.training)
             return predictor.get_post_fcov(self.training, X, diag)
         else:
             self._update_predictor()
             return self.predictor.get_post_fcov(self.training, X, diag)
 
+    def get_kernel_length_scale(self):
+        """
+        Return the Gaussian kernel length scale(s) (width) of the predictor.
+
+        With ARD, returns one length scale per input dimension; otherwise a single value.
+
+        Returns
+        -------
+        numpy.ndarray or None
+            Length scale(s). Shape (num_dim,) when ARD is used, (1,) otherwise.
+            None if the predictor is not set or not a GP with Gaussian kernel.
+        """
+        if self.predictor is None:
+            return None
+        self._update_predictor()
+        try:
+            cov = self.predictor.model.prior.cov
+        except AttributeError:
+            return None
+        if not hasattr(cov, "width"):
+            return None
+        return np.atleast_1d(np.asarray(cov.width).flatten())
+
+    def get_num_dim(self):
+        """
+        Return the input dimension (number of features) of the search space.
+
+        Returns
+        -------
+        int or None
+            The number of dimensions. None if neither training nor test data is set.
+        """
+        if (
+            self.training.X is not None
+            and self.training.X.shape[0] > 0
+        ):
+            return self.training.X.shape[1]
+        if self.test.X is not None and self.test.X.shape[0] > 0:
+            return self.test.X.shape[1]
+        return None
+
     def get_permutation_importance(self, n_perm: int, split_features_parallel=False):
         """
         Calculating permutation importance of model
@@ -508,7 +556,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False)
 
         if self.predictor is None:
             self._warn_no_predictor("get_post_fmean()")
-            predictor = gp_predictor(self.config)
+            predictor = self._make_gp_predictor()
             predictor.fit(self.training, 0)
             predictor.prepare(self.training)
             return predictor.get_permutation_importance(
@@ -588,7 +636,7 @@ def get_score(
         if predictor is None:
             if self.predictor is None:
                 self._warn_no_predictor("get_score()")
-                predictor = gp_predictor(self.config)
+                predictor = self._make_gp_predictor()
                 predictor.fit(training, 0, comm=self.mpicomm)
                 predictor.prepare(training)
             else:
@@ -899,9 +947,23 @@ def _init_predictor(self, is_rand_expans):
             If false, physbo.gp.Predictor is selected.
         """
         if is_rand_expans:
-            self.predictor = blm_predictor(self.config)
+            self.predictor = self._make_blm_predictor()
         else:
-            self.predictor = gp_predictor(self.config)
+            self.predictor = self._make_gp_predictor()
+
+    def _make_gp_predictor(self):
+        """Create a GP predictor, with ARD if self.ard is True."""
+        ard = self.ard
+        num_dim = self.get_num_dim()
+        model = gp.core.Model.create_default(ard=ard, num_dim=num_dim)
+        return gp_predictor(self.config, model=model)
+
+    def _make_blm_predictor(self):
+        """Create a BLM predictor, with ARD if self.ard is True."""
+        ard = self.ard
+        num_dim = self.get_num_dim()
+        model = gp.core.Model.create_default(ard=ard, num_dim=num_dim)
+        return blm_predictor(self.config, model=model)
 
     def _learn_hyperparameter(self, num_rand_basis):
         self.predictor.fit(self.training, num_rand_basis, comm=self.mpicomm)