Coverage for geometric_kernels / spaces / hypercube_graph.py: 94%
90 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-15 13:41 +0000
1"""
2This module provides the :class:`HypercubeGraph` space and the respective
3:class:`~.eigenfunctions.Eigenfunctions` subclass :class:`WalshFunctions`.
4"""
6from itertools import combinations
7from math import comb
9import lab as B
10import numpy as np
11from beartype.typing import List, Optional
13from geometric_kernels.lab_extras import dtype_double, float_like
14from geometric_kernels.spaces.base import DiscreteSpectrumSpace
15from geometric_kernels.spaces.eigenfunctions import (
16 Eigenfunctions,
17 EigenfunctionsWithAdditionTheorem,
18)
19from geometric_kernels.utils.special_functions import (
20 generalized_kravchuk_normalized,
21 walsh_function,
22)
23from geometric_kernels.utils.utils import chain, hamming_distance, log_binomial
26class WalshFunctions(EigenfunctionsWithAdditionTheorem):
27 r"""
28 Eigenfunctions of graph Laplacian on the hypercube graph $C^d$ whose nodes
29 are index by binary vectors in $\{0, 1\}^d$ are the Walsh
30 functions $w_T: C^d \to \{-1, 1\}$ given by
32 .. math:: w_T(x_0, .., x_{d-1}) = (-1)^{\sum_{i \in T} x_i},
34 enumerated by all possible subsets $T$ of the set $\{0, .., d-1\}$.
36 Levels are the whole eigenspaces, comprising all Walsh functions $w_T$ with
37 the same cardinality of $T$. The addition theorem for these is based on
38 certain discrete orthogonal polynomials called Kravchuk polynomials.
40 :param dim:
41 Dimension $d$ of the hypercube graph.
43 :param num_levels:
44 Specifies the number of levels of the Walsh functions.
45 """
47 def __init__(self, dim: int, num_levels: int) -> None:
48 if num_levels > dim + 1:
49 raise ValueError("The number of levels should be at most `dim`+1.")
50 self.dim = dim
51 self._num_levels = num_levels
52 self._num_eigenfunctions: Optional[int] = None # To be computed when needed.
54 def __call__(self, X: B.Bool, **kwargs) -> B.Float:
55 return B.stack(
56 *[
57 walsh_function(self.dim, list(cur_combination), X)
58 for level in range(self.num_levels)
59 for cur_combination in combinations(range(self.dim), level)
60 ],
61 axis=1,
62 )
64 def _addition_theorem(
65 self, X: B.Numeric, X2: Optional[B.Numeric] = None, **kwargs
66 ) -> B.Numeric:
68 if X2 is None:
69 X2 = X
71 hamming_distances = hamming_distance(X, X2)
73 values = []
75 kravchuk_normalized_j_minus_1, kravchuk_normalized_j_minus_2 = None, None
76 for level in range(self.num_levels):
77 cur_kravchuk_normalized = generalized_kravchuk_normalized(
78 self.dim,
79 level,
80 hamming_distances,
81 2,
82 kravchuk_normalized_j_minus_1,
83 kravchuk_normalized_j_minus_2,
84 ) # [N, N2]
85 kravchuk_normalized_j_minus_2 = kravchuk_normalized_j_minus_1
86 kravchuk_normalized_j_minus_1 = cur_kravchuk_normalized
88 values.append(
89 comb(self.dim, level) * cur_kravchuk_normalized[..., None] # [N, N2, 1]
90 )
92 return B.concat(*values, axis=-1) # [N, N2, L]
94 def _addition_theorem_diag(self, X: B.Numeric, **kwargs) -> B.Numeric:
95 """
96 These are certain easy to compute constants.
97 """
98 values = [
99 comb(self.dim, level) * B.ones(float_like(X), *X.shape[:-1], 1) # [N, 1]
100 for level in range(self.num_levels)
101 ]
102 return B.concat(*values, axis=1) # [N, L]
104 def weighted_outerproduct(
105 self,
106 weights: B.Numeric,
107 X: B.Numeric,
108 X2: Optional[B.Numeric] = None, # type: ignore
109 **kwargs,
110 ) -> B.Numeric:
111 if X2 is None:
112 X2 = X
114 hamming_distances = hamming_distance(X, X2)
116 result = B.zeros(B.dtype(weights), X.shape[0], X2.shape[0]) # [N, N2]
117 kravchuk_normalized_j_minus_1, kravchuk_normalized_j_minus_2 = None, None
118 for level in range(self.num_levels):
119 cur_kravchuk_normalized = generalized_kravchuk_normalized(
120 self.dim,
121 level,
122 hamming_distances,
123 2,
124 kravchuk_normalized_j_minus_1,
125 kravchuk_normalized_j_minus_2,
126 )
127 kravchuk_normalized_j_minus_2 = kravchuk_normalized_j_minus_1
128 kravchuk_normalized_j_minus_1 = cur_kravchuk_normalized
130 # Instead of multiplying weights by binomial coefficients, we sum their
131 # logs and then exponentiate the result for numerical stability.
132 # Furthermore, we save the computed Kravchuk polynomials for next iterations.
133 result += (
134 B.exp(B.log(weights[level]) + log_binomial(self.dim, level))
135 * cur_kravchuk_normalized
136 )
138 return result # [N, N2]
140 def weighted_outerproduct_diag(
141 self, weights: B.Numeric, X: B.Numeric, **kwargs
142 ) -> B.Numeric:
144 # Instead of multiplying weights by binomial coefficients, we sum their
145 # logs and then exponentiate the result for numerical stability.
146 result = sum(
147 B.exp(B.log(weights[level]) + log_binomial(self.dim, level))
148 * B.ones(float_like(X), *X.shape[:-1], 1)
149 for level in range(self.num_levels)
150 ) # [N, 1]
152 return B.reshape(result, *result.shape[:-1]) # [N,]
154 @property
155 def num_eigenfunctions(self) -> int:
156 if self._num_eigenfunctions is None:
157 self._num_eigenfunctions = sum(self.num_eigenfunctions_per_level)
158 return self._num_eigenfunctions
160 @property
161 def num_levels(self) -> int:
162 return self._num_levels
164 @property
165 def num_eigenfunctions_per_level(self) -> List[int]:
166 return [comb(self.dim, level) for level in range(self.num_levels)]
169class HypercubeGraph(DiscreteSpectrumSpace):
170 r"""
171 The GeometricKernels space representing the d-dimensional hypercube graph
172 $C^d = \{0, 1\}^d$, the combinatorial space of binary vectors of length $d$.
174 The elements of this space are represented by d-dimensional boolean vectors.
176 Levels are the whole eigenspaces.
178 .. note::
179 A tutorial on how to use this space is available in the
180 :doc:`HypercubeGraph.ipynb </examples/HypercubeGraph>` notebook.
182 .. note::
183 Since the degree matrix is a constant multiple of the identity, all
184 types of the graph Laplacian coincide on the hypercube graph up to a
185 constant, we choose the normalized Laplacian for numerical stability.
187 :param dim:
188 Dimension $d$ of the hypercube graph $C^d$, a positive integer.
190 .. admonition:: Citation
192 If you use this GeometricKernels space in your research, please consider
193 citing :cite:t:`borovitskiy2023`.
194 """
196 def __init__(self, dim: int):
197 if dim < 1:
198 raise ValueError("dim must be a positive integer.")
199 self.dim = dim
201 def __str__(self):
202 return f"HypercubeGraph({self.dim})"
204 @property
205 def dimension(self) -> int:
206 """
207 Returns d, the `dim` parameter that was passed down to `__init__`.
209 .. note:
210 Although this is a graph, and graphs are generally treated as
211 0-dimensional throughout GeometricKernels, we make an exception for
212 HypercubeGraph. This is because it helps maintain good behavior of
213 Matérn kernels with the usual values of the smoothness parameter
214 nu, i.e. nu = 1/2, nu = 3/2, nu = 5/2.
215 """
216 return self.dim
218 def get_eigenfunctions(self, num: int) -> Eigenfunctions:
219 """
220 Returns the :class:`~.WalshFunctions` object with `num` levels.
222 :param num:
223 Number of levels.
224 """
225 return WalshFunctions(self.dim, num)
227 def get_eigenvalues(self, num: int) -> B.Numeric:
228 eigenvalues = np.array(
229 [
230 2
231 * level
232 / self.dim # we assume normalized Laplacian (for numerical stability)
233 for level in range(num)
234 ]
235 )
236 return B.reshape(eigenvalues, -1, 1) # [num, 1]
238 def get_repeated_eigenvalues(self, num: int) -> B.Numeric:
239 eigenvalues_per_level = self.get_eigenvalues(num)
241 eigenfunctions = WalshFunctions(self.dim, num)
242 eigenvalues = chain(
243 B.squeeze(eigenvalues_per_level),
244 eigenfunctions.num_eigenfunctions_per_level,
245 ) # [J,]
246 return B.reshape(eigenvalues, -1, 1) # [J, 1]
248 def random(self, key: B.RandomState, number: int) -> B.Numeric:
249 r"""
250 Sample uniformly random points on the hypercube graph $C^d$.
252 Always returns [N, D] boolean array of the `key`'s backend.
254 :param key:
255 Either `np.random.RandomState`, `tf.random.Generator`,
256 `torch.Generator` or `jax.tensor` (representing random state).
257 :param number:
258 Number N of samples to draw.
260 :return:
261 An array of `number` uniformly random samples on the space.
262 """
263 key, random_points = B.random.rand(
264 key, dtype_double(key), number, self.dimension
265 )
267 random_points = random_points < 0.5
269 return key, random_points
271 @property
272 def element_shape(self):
273 """
274 :return:
275 [d].
276 """
277 return [self.dimension]
279 @property
280 def element_dtype(self):
281 """
282 :return:
283 B.Bool.
284 """
285 return B.Bool