Coverage for geometric_kernels/spaces/hypercube_graph.py: 94%
90 statements
« prev ^ index » next coverage.py v7.11.3, created at 2025-11-16 21:43 +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 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 = kravchuk_normalized(
78 self.dim,
79 level,
80 hamming_distances,
81 kravchuk_normalized_j_minus_1,
82 kravchuk_normalized_j_minus_2,
83 ) # [N, N2]
84 kravchuk_normalized_j_minus_2 = kravchuk_normalized_j_minus_1
85 kravchuk_normalized_j_minus_1 = cur_kravchuk_normalized
87 values.append(
88 comb(self.dim, level) * cur_kravchuk_normalized[..., None] # [N, N2, 1]
89 )
91 return B.concat(*values, axis=-1) # [N, N2, L]
93 def _addition_theorem_diag(self, X: B.Numeric, **kwargs) -> B.Numeric:
94 """
95 These are certain easy to compute constants.
96 """
97 values = [
98 comb(self.dim, level) * B.ones(float_like(X), *X.shape[:-1], 1) # [N, 1]
99 for level in range(self.num_levels)
100 ]
101 return B.concat(*values, axis=1) # [N, L]
103 def weighted_outerproduct(
104 self,
105 weights: B.Numeric,
106 X: B.Numeric,
107 X2: Optional[B.Numeric] = None, # type: ignore
108 **kwargs,
109 ) -> B.Numeric:
110 if X2 is None:
111 X2 = X
113 hamming_distances = hamming_distance(X, X2)
115 result = B.zeros(B.dtype(weights), X.shape[0], X2.shape[0]) # [N, N2]
116 kravchuk_normalized_j_minus_1, kravchuk_normalized_j_minus_2 = None, None
117 for level in range(self.num_levels):
118 cur_kravchuk_normalized = kravchuk_normalized(
119 self.dim,
120 level,
121 hamming_distances,
122 kravchuk_normalized_j_minus_1,
123 kravchuk_normalized_j_minus_2,
124 )
125 kravchuk_normalized_j_minus_2 = kravchuk_normalized_j_minus_1
126 kravchuk_normalized_j_minus_1 = cur_kravchuk_normalized
128 # Instead of multiplying weights by binomial coefficients, we sum their
129 # logs and then exponentiate the result for numerical stability.
130 # Furthermore, we save the computed Kravchuk polynomials for next iterations.
131 result += (
132 B.exp(B.log(weights[level]) + log_binomial(self.dim, level))
133 * cur_kravchuk_normalized
134 )
136 return result # [N, N2]
138 def weighted_outerproduct_diag(
139 self, weights: B.Numeric, X: B.Numeric, **kwargs
140 ) -> B.Numeric:
142 # Instead of multiplying weights by binomial coefficients, we sum their
143 # logs and then exponentiate the result for numerical stability.
144 result = sum(
145 B.exp(B.log(weights[level]) + log_binomial(self.dim, level))
146 * B.ones(float_like(X), *X.shape[:-1], 1)
147 for level in range(self.num_levels)
148 ) # [N, 1]
150 return B.reshape(result, *result.shape[:-1]) # [N,]
152 @property
153 def num_eigenfunctions(self) -> int:
154 if self._num_eigenfunctions is None:
155 self._num_eigenfunctions = sum(self.num_eigenfunctions_per_level)
156 return self._num_eigenfunctions
158 @property
159 def num_levels(self) -> int:
160 return self._num_levels
162 @property
163 def num_eigenfunctions_per_level(self) -> List[int]:
164 return [comb(self.dim, level) for level in range(self.num_levels)]
167class HypercubeGraph(DiscreteSpectrumSpace):
168 r"""
169 The GeometricKernels space representing the d-dimensional hypercube graph
170 $C^d = \{0, 1\}^d$, the combinatorial space of binary vectors of length $d$.
172 The elements of this space are represented by d-dimensional boolean vectors.
174 Levels are the whole eigenspaces.
176 .. note::
177 A tutorial on how to use this space is available in the
178 :doc:`HypercubeGraph.ipynb </examples/HypercubeGraph>` notebook.
180 .. note::
181 Since the degree matrix is a constant multiple of the identity, all
182 types of the graph Laplacian coincide on the hypercube graph up to a
183 constant, we choose the normalized Laplacian for numerical stability.
185 :param dim:
186 Dimension $d$ of the hypercube graph $C^d$, a positive integer.
188 .. admonition:: Citation
190 If you use this GeometricKernels space in your research, please consider
191 citing :cite:t:`borovitskiy2023`.
192 """
194 def __init__(self, dim: int):
195 if dim < 1:
196 raise ValueError("dim must be a positive integer.")
197 self.dim = dim
199 def __str__(self):
200 return f"HypercubeGraph({self.dim})"
202 @property
203 def dimension(self) -> int:
204 """
205 Returns d, the `dim` parameter that was passed down to `__init__`.
207 .. note:
208 Although this is a graph, and graphs are generally treated as
209 0-dimensional throughout GeometricKernels, we make an exception for
210 HypercubeGraph. This is because it helps maintain good behavior of
211 Matérn kernels with the usual values of the smoothness parameter
212 nu, i.e. nu = 1/2, nu = 3/2, nu = 5/2.
213 """
214 return self.dim
216 def get_eigenfunctions(self, num: int) -> Eigenfunctions:
217 """
218 Returns the :class:`~.WalshFunctions` object with `num` levels.
220 :param num:
221 Number of levels.
222 """
223 return WalshFunctions(self.dim, num)
225 def get_eigenvalues(self, num: int) -> B.Numeric:
226 eigenvalues = np.array(
227 [
228 2
229 * level
230 / self.dim # we assume normalized Laplacian (for numerical stability)
231 for level in range(num)
232 ]
233 )
234 return B.reshape(eigenvalues, -1, 1) # [num, 1]
236 def get_repeated_eigenvalues(self, num: int) -> B.Numeric:
237 eigenvalues_per_level = self.get_eigenvalues(num)
239 eigenfunctions = WalshFunctions(self.dim, num)
240 eigenvalues = chain(
241 B.squeeze(eigenvalues_per_level),
242 eigenfunctions.num_eigenfunctions_per_level,
243 ) # [J,]
244 return B.reshape(eigenvalues, -1, 1) # [J, 1]
246 def random(self, key: B.RandomState, number: int) -> B.Numeric:
247 r"""
248 Sample uniformly random points on the hypercube graph $C^d$.
250 Always returns [N, D] boolean array of the `key`'s backend.
252 :param key:
253 Either `np.random.RandomState`, `tf.random.Generator`,
254 `torch.Generator` or `jax.tensor` (representing random state).
255 :param number:
256 Number N of samples to draw.
258 :return:
259 An array of `number` uniformly random samples on the space.
260 """
261 key, random_points = B.random.rand(
262 key, dtype_double(key), number, self.dimension
263 )
265 random_points = random_points < 0.5
267 return key, random_points
269 @property
270 def element_shape(self):
271 """
272 :return:
273 [d].
274 """
275 return [self.dimension]
277 @property
278 def element_dtype(self):
279 """
280 :return:
281 B.Bool.
282 """
283 return B.Bool