Hypergraph Lumpability and Bayesian Refinement

Download original file

This is the original source text. Reading it does not run or verify the program.

Read the typeset paper

1 """Finite reference implementation for interaction quotients resolved by projectors.
2 
3 Tensor convention:
4 A multilinear map of order m, A : R^n x ... x R^n -> R^p, is represented by
5 an array with shape ``(p, n, ..., n)``.  Axis 0 is the output axis and axes
6 1..m are the ordered input slots.
7 
8 The implementation is intentionally dense and explicit.  It is meant for
9 small certificates attached to theorems, not production workloads at large scale.
10 """
11 
12 from __future__ import annotations
13 
14 from dataclasses import dataclass
15 from itertools import product
16 from math import exp
17 from typing import Callable, Iterable, Sequence
18 
19 import numpy as np
20 
21 Array = np.ndarray
22 
23 
24 class PRIQError(ValueError):
25     """Raised when a finite PRIQ object violates its required shape contract."""
26 
27 
28 def _as_float_array(value: Array | Sequence[float], *, ndim: int | None = None) -> Array:
29     arr = np.asarray(value, dtype=float)
30     if ndim is not None and arr.ndim != ndim:
31         raise PRIQError(f"expected an array with ndim={ndim}, received shape {arr.shape}")
32     if not np.all(np.isfinite(arr)):
33         raise PRIQError("all entries must be finite")
34     return arr
35 
36 
37 def block_sum_matrix(packet_map: Sequence[int], coarse_count: int | None = None) -> Array:
38     """Return the block sum aggregation ``C`` determined by a packet map.
39 
40     ``packet_map[j]`` is the coarse packet containing fine coordinate ``j``.
41     The returned matrix satisfies ``(C @ x)[a] = sum_{j:pi(j)=a} x[j]``.
42     """
43 
44     pi = np.asarray(packet_map, dtype=int)
45     if pi.ndim != 1 or pi.size == 0:
46         raise PRIQError("packet_map must be a nonempty sequence of one dimension")
47     if np.any(pi < 0):
48         raise PRIQError("packet labels must be nonnegative")
49     q = int(pi.max()) + 1 if coarse_count is None else int(coarse_count)
50     if q <= int(pi.max()):
51         raise PRIQError("coarse_count does not contain every packet label")
52     C = np.zeros((q, pi.size), dtype=float)
53     C[pi, np.arange(pi.size)] = 1.0
54     if np.any(C.sum(axis=1) == 0):
55         raise PRIQError("every coarse packet must have a nonempty fine fiber")
56     return C
57 
58 
59 def averaging_section(packet_map: Sequence[int], coarse_count: int | None = None) -> Array:
60     """Return the canonical right inverse ``E`` with equal weights in each block."""
61 
62     pi = np.asarray(packet_map, dtype=int)
63     C = block_sum_matrix(pi, coarse_count)
64     q, n = C.shape
65     E = np.zeros((n, q), dtype=float)
66     for a in range(q):
67         members = np.flatnonzero(pi == a)
68         E[members, a] = 1.0 / members.size
69     if not np.allclose(C @ E, np.eye(q), atol=1e-14, rtol=1e-14):
70         raise PRIQError("failed to construct a right inverse")
71     return E
72 
73 
74 def kernel_basis(C: Array, *, rtol: float = 1e-12) -> Array:
75     """Compute an orthonormal basis for ``ker(C)`` using an SVD."""
76 
77     C = _as_float_array(C, ndim=2)
78     _, singular, vh = np.linalg.svd(C, full_matrices=True)
79     if singular.size == 0:
80         return np.eye(C.shape[1])
81     threshold = rtol * max(C.shape) * singular.max(initial=0.0)
82     rank = int(np.sum(singular > threshold))
83     return vh[rank:].T.copy()
84 
85 
86 def validate_tensor(tensor: Array) -> tuple[int, int, int]:
87     """Validate a multilinear tensor whose input dimensions are equal.
88 
89     Returns ``(order, output_dimension, input_dimension)``.
90     """
91 
92     A = _as_float_array(tensor)
93     if A.ndim < 2:
94         raise PRIQError("a multilinear tensor must have one output and at least one input axis")
95     m = A.ndim - 1
96     input_dims = A.shape[1:]
97     if len(set(input_dims)) != 1:
98         raise PRIQError(f"all input dimensions must agree, received {input_dims}")
99     return m, A.shape[0], input_dims[0]
100 
101 
102 def transform_tensor(
103     tensor: Array,
104     output_map: Array,
105     input_maps: Sequence[Array],
106 ) -> Array:
107     """Apply a linear map to the output and one linear map to every input slot.
108 
109     If ``tensor`` has coefficients ``A[o,i1,...,im]``, ``output_map`` has shape
110     ``(O,o)``, and input map ``r`` has shape ``(ir,Ir)``, the returned tensor is
111 
112     ``B[O,I1,...,Im] = output_map[O,o] A[o,i1,...,im] Π_r input_map_r[ir,Ir]``.
113     """
114 
115     A = _as_float_array(tensor)
116     if A.ndim < 2:
117         raise PRIQError("tensor must have at least one input slot")
118     m = A.ndim - 1
119     if len(input_maps) != m:
120         raise PRIQError(f"expected {m} input maps, received {len(input_maps)}")
121     O = _as_float_array(output_map, ndim=2)
122     if O.shape[1] != A.shape[0]:
123         raise PRIQError("output_map has an incompatible source output dimension")
124 
125     # NumPy's integer subscript form of einsum avoids fragile generated strings.
126     # It supports every modest tensor order used by this package.
127     old_output = 0
128     old_inputs = list(range(1, m + 1))
129     new_output = m + 1
130     new_inputs = list(range(m + 2, 2 * m + 2))
131 
132     operands: list[object] = [A, [old_output, *old_inputs], O, [new_output, old_output]]
133     for axis, (M, old_label, new_label) in enumerate(
134         zip(input_maps, old_inputs, new_inputs, strict=True), start=1
135     ):
136         M_arr = _as_float_array(M, ndim=2)
137         if M_arr.shape[0] != A.shape[axis]:
138             raise PRIQError(f"input map for slot {axis} has an incompatible source input dimension")
139         operands.extend([M_arr, [old_label, new_label]])
140     operands.append([new_output, *new_inputs])
141     return np.einsum(*operands, optimize=True)
142 
143 
144 def apply_tensor(tensor: Array, *vectors: Array) -> Array:
145     """Evaluate a multilinear tensor on ordered input vectors."""
146 
147     A = _as_float_array(tensor)
148     m, _, input_dim = validate_tensor(A)
149     if len(vectors) != m:
150         raise PRIQError(f"expected {m} vectors, received {len(vectors)}")
151     operands: list[object] = [A, list(range(m + 1))]
152     for slot, vector in enumerate(vectors, start=1):
153         x = _as_float_array(vector, ndim=1)
154         if x.size != input_dim:
155             raise PRIQError(f"vector in slot {slot} has length {x.size}, expected {input_dim}")
156         operands.extend([x, [slot]])
157     operands.append([0])
158     return np.einsum(*operands, optimize=True)
159 
160 
161 def symmetrize_inputs(tensor: Array) -> Array:
162     """Average a tensor over all permutations of its input axes."""
163 
164     from itertools import permutations
165 
166     A = _as_float_array(tensor)
167     m, _, _ = validate_tensor(A)
168     result = np.zeros_like(A)
169     perms = list(permutations(range(1, m + 1)))
170     for perm in perms:
171         result += np.transpose(A, axes=(0, *perm))
172     return result / len(perms)
173 
174 
175 def quotient_defect(fine_tensor: Array, coarse_tensor: Array, C: Array) -> Array:
176     """Return ``C A_f - A_c C^{x m}`` as a coefficient tensor."""
177 
178     Af = _as_float_array(fine_tensor)
179     Ac = _as_float_array(coarse_tensor)
180     C = _as_float_array(C, ndim=2)
181     m_f, out_f, in_f = validate_tensor(Af)
182     m_c, out_c, in_c = validate_tensor(Ac)
183     if m_f != m_c:
184         raise PRIQError("fine and coarse tensors must have the same order")
185     if C.shape != (out_c, out_f) or in_f != out_f or in_c != out_c:
186         raise PRIQError("this reference routine expects square fine/coarse state tensors")
187     left = transform_tensor(Af, C, [np.eye(in_f)] * m_f)
188     right = transform_tensor(Ac, np.eye(out_c), [C] * m_f)
189     return left - right
190 
191 
192 def exact_lift_tensor(
193     coarse_tensor: Array,
194     C: Array,
195     E: Array,
196     *,
197     hidden_tensor: Array | None = None,
198 ) -> Array:
199     """Lift a coarse tensor and optionally add a tensor with hidden output.
200 
201     The base lift is ``E A_c C^{x m}``.  Any supplied hidden tensor is projected
202     through ``Q = I - E C`` on its output, so its aggregate is identically zero.
203     """
204 
205     Ac = _as_float_array(coarse_tensor)
206     C = _as_float_array(C, ndim=2)
207     E = _as_float_array(E, ndim=2)
208     m, q_out, q_in = validate_tensor(Ac)
209     q, n = C.shape
210     if q_out != q or q_in != q or E.shape != (n, q):
211         raise PRIQError("coarse tensor, aggregation, and section have incompatible dimensions")
212     if not np.allclose(C @ E, np.eye(q), atol=1e-10, rtol=1e-10):
213         raise PRIQError("E must be a right inverse of C")
214     base = transform_tensor(Ac, E, [C] * m)
215     if hidden_tensor is None:
216         return base
217     H = _as_float_array(hidden_tensor)
218     if H.shape != base.shape:
219         raise PRIQError(f"hidden tensor has shape {H.shape}, expected {base.shape}")
220     Q = np.eye(n) - E @ C
221     hidden = transform_tensor(H, Q, [np.eye(n)] * m)
222     return base + hidden
223 
224 
225 def aggregate_output_coefficients(fine_tensor: Array, C: Array) -> Array:
226     """Return the coefficient array ``C A_f`` after output aggregation."""
227 
228     Af = _as_float_array(fine_tensor)
229     C = _as_float_array(C, ndim=2)
230     m, out_dim, in_dim = validate_tensor(Af)
231     if C.shape[1] != out_dim:
232         raise PRIQError("aggregation has incompatible output dimension")
233     return transform_tensor(Af, C, [np.eye(in_dim)] * m)
234 
235 
236 @dataclass(frozen=True)
237 class RecoveryResult:
238     coarse_tensor: Array
239     max_fiber_variation: float
240     witness: tuple[int, tuple[int, ...], tuple[int, ...]] | None
241 
242 
243 def recover_coarse_tensor(
244     fine_tensor: Array,
245     packet_map: Sequence[int],
246     *,
247     tolerance: float = 1e-10,
248     representative: str = "mean",
249 ) -> RecoveryResult:
250     """Recover a coarse tensor from fiber values in packet coordinates.
251 
252     ``representative`` may be ``"mean"`` or ``"midrange"``.  The returned
253     variation records the full spread within the worst fiber.  A numerical run
254     may accept a spread below ``tolerance``.  Algebraic exactness still means
255     that the spread is zero.
256     """
257 
258     Af = _as_float_array(fine_tensor)
259     m, n_out, n_in = validate_tensor(Af)
260     if n_out != n_in:
261         raise PRIQError("block sum recovery expects a state tensor with equal input and output dimensions")
262     pi = np.asarray(packet_map, dtype=int)
263     if pi.size != n_in:
264         raise PRIQError("packet_map length must equal the fine state dimension")
265     C = block_sum_matrix(pi)
266     q = C.shape[0]
267     G = aggregate_output_coefficients(Af, C)
268     fibers = [np.flatnonzero(pi == a).tolist() for a in range(q)]
269     coarse = np.empty((q, *([q] * m)), dtype=float)
270     max_variation = 0.0
271     witness: tuple[int, tuple[int, ...], tuple[int, ...]] | None = None
272 
273     for coarse_inputs in product(range(q), repeat=m):
274         fine_tuples = list(product(*(fibers[a] for a in coarse_inputs)))
275         for a in range(q):
276             values = np.asarray([G[(a, *idx)] for idx in fine_tuples], dtype=float)
277             variation = float(values.max() - values.min())
278             if variation > max_variation:
279                 max_variation = variation
280                 i_min = fine_tuples[int(np.argmin(values))]
281                 i_max = fine_tuples[int(np.argmax(values))]
282                 witness = (a, tuple(i_min), tuple(i_max))
283             if representative == "mean":
284                 value = float(values.mean())
285             elif representative == "midrange":
286                 value = float(0.5 * (values.max() + values.min()))
287             else:
288                 raise PRIQError("representative must be 'mean' or 'midrange'")
289             coarse[(a, *coarse_inputs)] = value
290 
291     if max_variation <= tolerance:
292         witness = None
293     return RecoveryResult(coarse, max_variation, witness)
294 
295 
296 def slot_kernel_residual(fine_tensor: Array, C: Array) -> float:
297     """Return the largest Frobenius residual after inserting ``ker(C)`` in one slot."""
298 
299     Af = _as_float_array(fine_tensor)
300     C = _as_float_array(C, ndim=2)
301     m, n_out, n_in = validate_tensor(Af)
302     if C.shape[1] != n_out or n_in != n_out:
303         raise PRIQError("the slot kernel test expects a square fine state tensor")
304     K = kernel_basis(C)
305     if K.shape[1] == 0:
306         return 0.0
307     residual = 0.0
308     for slot in range(m):
309         maps = [np.eye(n_in) for _ in range(m)]
310         maps[slot] = K
311         transformed = transform_tensor(Af, C, maps)
312         residual = max(residual, float(np.linalg.norm(transformed)))
313     return residual
314 
315 
316 def bayes_update(mu: Array, likelihood: Array) -> Array:
317     """Return the normalized positive likelihood update."""
318 
319     p = _as_float_array(mu, ndim=1)
320     lam = _as_float_array(likelihood, ndim=1)
321     if p.size != lam.size:
322         raise PRIQError("state and likelihood dimensions differ")
323     if np.any(p < 0) or np.any(lam < 0):
324         raise PRIQError("Bayesian update requires nonnegative state and likelihood")
325     unnormalized = lam * p
326     evidence = float(unnormalized.sum())
327     if evidence <= 0:
328         raise PRIQError("Bayesian evidence must be positive")
329     return unnormalized / evidence
330 
331 
332 def bayes_commutation_residual(
333     mu_fine: Array,
334     lambda_fine: Array,
335     lambda_coarse: Array,
336     C: Array,
337     *,
338     ord: int | float = 1,
339 ) -> float:
340     """Compute ``|| C B_f(mu) - B_c(C mu) ||``."""
341 
342     C = _as_float_array(C, ndim=2)
343     lhs = C @ bayes_update(mu_fine, lambda_fine)
344     rhs = bayes_update(C @ _as_float_array(mu_fine, ndim=1), lambda_coarse)
345     return float(np.linalg.norm(lhs - rhs, ord=ord))
346 
347 
348 def normalization_perturbation_bound(u: Array, v: Array) -> float:
349     """Return the symmetric l1 bound for normalization of positive vectors."""
350 
351     u = _as_float_array(u, ndim=1)
352     v = _as_float_array(v, ndim=1)
353     if np.any(u < 0) or np.any(v < 0):
354         raise PRIQError("normalization vectors must be nonnegative")
355     a = float(np.sum(u))
356     b = float(np.sum(v))
357     if a <= 0 or b <= 0:
358         raise PRIQError("normalization masses must be positive")
359     return 2.0 * float(np.linalg.norm(u - v, ord=1)) / max(a, b)
360 
361 
362 def rk4_step(field: Callable[[Array], Array], state: Array, step: float) -> Array:
363     """Take one deterministic classical RK4 step."""
364 
365     x = _as_float_array(state, ndim=1)
366     h = float(step)
367     if h <= 0:
368         raise PRIQError("step must be positive")
369     k1 = _as_float_array(field(x), ndim=1)
370     k2 = _as_float_array(field(x + 0.5 * h * k1), ndim=1)
371     k3 = _as_float_array(field(x + 0.5 * h * k2), ndim=1)
372     k4 = _as_float_array(field(x + h * k3), ndim=1)
373     if not (k1.size == k2.size == k3.size == k4.size == x.size):
374         raise PRIQError("field returned an incompatible state dimension")
375     return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
376 
377 
378 def integrate_rk4(
379     field: Callable[[Array], Array],
380     initial: Array,
381     duration: float,
382     *,
383     max_step: float = 1e-3,
384     project_simplex: bool = False,
385 ) -> Array:
386     """Integrate an autonomous finite ODE over one interval."""
387 
388     T = float(duration)
389     if T < 0:
390         raise PRIQError("duration must be nonnegative")
391     x = _as_float_array(initial, ndim=1).copy()
392     if T == 0:
393         return x
394     steps = max(1, int(np.ceil(T / max_step)))
395     h = T / steps
396     for _ in range(steps):
397         x = rk4_step(field, x, h)
398         if project_simplex:
399             # Exact commutation tests do not use projection.  It is available only
400             # as a numerical safeguard for explicitly approximate runs.
401             x = np.maximum(x, 0.0)
402             total = float(x.sum())
403             if total <= 0:
404                 raise PRIQError("simplex projection encountered zero total mass")
405             x /= total
406     return x
407 
408 
409 @dataclass(frozen=True)
410 class FilterTrace:
411     fine_posteriors: Array
412     coarse_posteriors: Array
413     fine_predictions: Array
414     coarse_predictions: Array
415     fine_evidence: Array
416     coarse_evidence: Array
417 
418 
419 def simulate_prediction_update_filter(
420     fine_initial: Array,
421     coarse_initial: Array,
422     fine_field: Callable[[Array], Array],
423     coarse_field: Callable[[Array], Array],
424     fine_likelihoods: Sequence[Array],
425     coarse_likelihoods: Sequence[Array],
426     durations: Sequence[float],
427     *,
428     max_step: float = 1e-3,
429 ) -> FilterTrace:
430     """Run matched prediction and update recursions for the fine and coarse models."""
431 
432     if not (len(fine_likelihoods) == len(coarse_likelihoods) == len(durations)):
433         raise PRIQError("likelihood and duration sequences must have the same length")
434     xf = _as_float_array(fine_initial, ndim=1).copy()
435     xc = _as_float_array(coarse_initial, ndim=1).copy()
436     fine_posts = [xf.copy()]
437     coarse_posts = [xc.copy()]
438     fine_preds: list[Array] = []
439     coarse_preds: list[Array] = []
440     fine_evidence: list[float] = []
441     coarse_evidence: list[float] = []
442 
443     for lf, lc, T in zip(fine_likelihoods, coarse_likelihoods, durations, strict=True):
444         xf_hat = integrate_rk4(fine_field, xf, T, max_step=max_step)
445         xc_hat = integrate_rk4(coarse_field, xc, T, max_step=max_step)
446         lf_arr = _as_float_array(lf, ndim=1)
447         lc_arr = _as_float_array(lc, ndim=1)
448         zf = float(lf_arr @ xf_hat)
449         zc = float(lc_arr @ xc_hat)
450         xf = bayes_update(xf_hat, lf_arr)
451         xc = bayes_update(xc_hat, lc_arr)
452         fine_preds.append(xf_hat)
453         coarse_preds.append(xc_hat)
454         fine_evidence.append(zf)
455         coarse_evidence.append(zc)
456         fine_posts.append(xf.copy())
457         coarse_posts.append(xc.copy())
458 
459     return FilterTrace(
460         fine_posteriors=np.asarray(fine_posts),
461         coarse_posteriors=np.asarray(coarse_posts),
462         fine_predictions=np.asarray(fine_preds),
463         coarse_predictions=np.asarray(coarse_preds),
464         fine_evidence=np.asarray(fine_evidence),
465         coarse_evidence=np.asarray(coarse_evidence),
466     )
467 
468 
469 def posterior_error_certificate(
470     initial_error: float,
471     flow_lipschitz: Sequence[float],
472     bayes_lipschitz: Sequence[float],
473     flow_defects: Sequence[float],
474     bayes_defects: Sequence[float],
475 ) -> tuple[float, Array, Array]:
476     """Propagate the affine error recurrence for prediction and update."""
477 
478     arrays = [
479         np.asarray(flow_lipschitz, dtype=float),
480         np.asarray(bayes_lipschitz, dtype=float),
481         np.asarray(flow_defects, dtype=float),
482         np.asarray(bayes_defects, dtype=float),
483     ]
484     lengths = {arr.size for arr in arrays}
485     if len(lengths) != 1:
486         raise PRIQError("all certificate sequences must have the same length")
487     if any(np.any(arr < 0) for arr in arrays):
488         raise PRIQError("certificate factors and defects must be nonnegative")
489     e_plus = float(initial_error)
490     minus_trace: list[float] = []
491     plus_trace: list[float] = [e_plus]
492     for Lphi, LB, dphi, dB in zip(*arrays, strict=True):
493         e_minus = float(Lphi * e_plus + dphi)
494         e_plus = float(LB * e_minus + dB)
495         minus_trace.append(e_minus)
496         plus_trace.append(e_plus)
497     return e_plus, np.asarray(minus_trace), np.asarray(plus_trace)
498 
499 
500 def hidden_exchange_fields(
501     *,
502     a: float = 0.7,
503     b: float = 0.4,
504     r0: float = 0.25,
505     r1: float = 0.30,
506     s0: float = 0.20,
507     s1: float = 0.18,
508     u0: float = 0.22,
509     u1: float = 0.24,
510     v0: float = 0.19,
511     v1: float = 0.20,
512     epsilon: float = 0.0,
513 ) -> tuple[Callable[[Array], Array], Callable[[Array], Array], Array]:
514     """Return the exact example with four fine states and two coarse states."""
515 
516     params = [a, b, r0, r1, s0, s1, u0, u1, v0, v1, epsilon]
517     if any(value < 0 for value in params):
518         raise PRIQError("all example rates must be nonnegative")
519     C = block_sum_matrix([0, 0, 1, 1])
520 
521     def coarse(z: Array) -> Array:
522         z = _as_float_array(z, ndim=1)
523         return np.asarray([-a * z[0] + b * z[1], a * z[0] - b * z[1]])
524 
525     def fine(x: Array) -> Array:
526         x = _as_float_array(x, ndim=1)
527         h12 = (r0 + r1 * x[2]) * x[1] - (s0 + s1 * x[3]) * x[0]
528         h34 = (u0 + u1 * x[0]) * x[3] - (v0 + v1 * x[1]) * x[2]
529         dx = np.asarray(
530             [
531                 -a * x[0] + b * x[2] + h12,
532                 -a * x[1] + b * x[3] - h12,
533                 a * x[0] - b * x[2] + h34,
534                 a * x[1] - b * x[3] - h34,
535             ],
536             dtype=float,
537         )
538         if epsilon:
539             dx += epsilon * x[0] * np.asarray([-1.0, 0.0, 1.0, 0.0])
540         return dx
541 
542     return fine, coarse, C
543 
544 
545 def hyperedge_energy_statistics(coefficients: Array) -> tuple[Array, float, float]:
546     """Return normalized squared energies, Shannon entropy, and effective count."""
547 
548     values = _as_float_array(coefficients).ravel()
549     energy = values * values
550     total = float(energy.sum())
551     if total <= 0:
552         raise PRIQError("at least one hyperedge coefficient must be nonzero")
553     probabilities = energy / total
554     positive = probabilities > 0
555     entropy = float(-np.sum(probabilities[positive] * np.log(probabilities[positive])))
556     return probabilities, entropy, exp(entropy)
557 
558 
559 def best_k_residuals(coefficients: Array) -> Array:
560     """Return exact Frobenius residual after retaining each number of largest blocks."""
561 
562     values = _as_float_array(coefficients).ravel()
563     energies = np.sort(values * values)[::-1]
564     residual_sq = np.zeros(energies.size + 1, dtype=float)
565     residual_sq[:-1] = np.cumsum(energies[::-1])[::-1]
566     return np.sqrt(residual_sq)
567 
568 
569 def total_variation(p: Array, q: Array) -> float:
570     """Finite total variation distance."""
571 
572     p_arr = _as_float_array(p, ndim=1)
573     q_arr = _as_float_array(q, ndim=1)
574     if p_arr.size != q_arr.size:
575         raise PRIQError("probability vectors must have the same dimension")
576     return 0.5 * float(np.linalg.norm(p_arr - q_arr, ord=1))

Back to source controls

Back to Hypergraph Lumpability and Bayesian Refinement