Universal Coupling Prior

1. Mathematical Formalism

Bayesian Coupling Parametrisation

The SCPN coupling matrix \(K_{ij}\) is parametrised by two scalars:

\[K_{ij} = K_{\text{base}} \cdot e^{-\alpha_d \cdot |i - j|}\]

where: - \(K_{\text{base}}\) is the base coupling strength (nearest-neighbour) - \(\alpha_d\) is the distance-decay exponent - \(|i - j|\) is the topological distance between oscillators \(i\) and \(j\)

This exponential-decay model captures the empirical observation that coupling strength decreases with distance in most physical, biological, and engineered networks.

The Prior Distribution

The UniversalPrior encodes the empirical distribution of coupling parameters across 25 domain-specific configurations (R4-A3 cross-domain transfer study):

\[K_{\text{base}} \sim \mathcal{N}(\mu_K = 0.47, \sigma_K = 0.09)\]

\[\alpha_d \sim \mathcal{N}(\mu_\alpha = 0.25, \sigma_\alpha = 0.07)\]

The two parameters are assumed independent (diagonal covariance).

Unnormalised Log-Probability

The log_probability method computes:

\[\log p(K_{\text{base}}, \alpha_d) = -\frac{1}{2}\left(\frac{K_{\text{base}} - \mu_K}{\sigma_K}\right)^2 - \frac{1}{2}\left(\frac{\alpha_d - \mu_\alpha}{\sigma_\alpha}\right)^2\]

This is the unnormalised log-density of the bivariate Gaussian prior. The normalisation constant \(-\log(2\pi\sigma_K\sigma_\alpha)\) is omitted because it is constant across parameter evaluations and cancels in Bayesian posterior ratios.

Dörfler-Bullo Critical Coupling

The estimate_Kc method combines the prior with the critical coupling estimate from algebraic graph theory:

\[K_c = \frac{\max_{i,j} |\omega_i - \omega_j|}{\lambda_2(L(W))}\]

where \(L(W)\) is the graph Laplacian of the coupling matrix built from the prior's MAP estimate, and \(\lambda_2\) is its algebraic connectivity (Fiedler value).

This estimate gives the minimum coupling strength needed for frequency synchronisation given the natural frequency spread \(\max |\omega_i - \omega_j|\) and the network topology. The public boundary validates omegas as a finite one-dimensional real-valued frequency vector before constructing the distance-decay matrix or dispatching to the spectral critical-coupling calculation. Boolean aliases, complex values, non-finite values, empty vectors, and length mismatches are rejected. The alias boundary includes Python bool, NumPy boolean scalars, and object arrays containing either form so frequencies cannot be silently coerced into 0.0/1.0 weights.

Distance-Decay Matrix

The coupling matrix generated by the prior is:

\[W_{ij} = \begin{cases} K_{\text{base}} \cdot e^{-\alpha_d \cdot |i-j|} & i \neq j \\ 0 & i = j \end{cases}\]

Properties: - Symmetric: \(W_{ij} = W_{ji}\) (undirected coupling) - Non-negative: \(W_{ij} \geq 0\) for all \(i, j\) - Monotone decreasing: \(W_{ij} > W_{ik}\) when \(|i-j| < |i-k|\) - Zero diagonal: no self-coupling

Dimensionality Reduction

The prior collapses the coupling parametrisation from \(O(N^2)\) free parameters (full matrix) to 2 scalar parameters \((K_{\text{base}}, \alpha_d)\). Combined with the Dörfler-Bullo \(K_c\) constraint, the auto-tune search space is effectively:

\[\text{dim}(\text{search}) = 2 \quad \text{instead of} \quad N^2\]

For \(N = 256\), this reduces from 65,536 to 2 parameters.

---

2. Theoretical Context

Why a Universal Prior?

Different domains (EEG, power grids, plasma, chemical oscillators) have different optimal coupling parameters. However, analysing 25 domain-specific configurations revealed that the parameters cluster tightly around \(K_{\text{base}} \approx 0.47\) and \(\alpha_d \approx 0.25\).

This universality arises because: 1. Scale invariance: The normalised coupling \(K/K_c\) has similar optimal values across domains (typically \(K/K_c \in [1.5, 3.0]\)). 2. Small-world topology: Distance-decay with \(\alpha_d \approx 0.25\) produces small-world-like connectivity: strong local coupling with weak but non-zero long-range connections. 3. Stability margins: These parameter values place the system well above the synchronisation threshold but below the regime where numerical stiffness becomes problematic.

The 25-Domainpack Study

The empirical prior was derived from the R4-A3 cross-domain transfer analysis (Stankovski 2017). The 25 domains include:

• Neuroscience: EEG (alpha, beta, gamma bands), fMRI resting-state, MEG source-space
• Physics: coupled pendula, Josephson junction arrays, laser arrays
• Engineering: power grid frequency regulation, clock synchronisation
• Biology: circadian rhythms, cardiac pacemaker cells, firefly synchronisation
• Chemistry: Belousov-Zhabotinsky reaction, electrochemical oscillators

For each domain, the optimal \((K_{\text{base}}, \alpha_d)\) was determined by maximising the domain-specific objective (e.g., order parameter \(R\) for neuroscience, frequency deviation for power grids). The resulting distribution was well-approximated by the bivariate Gaussian.

Historical Context

• Stankovski, T. et al. (2017): "Coupling functions: Universal insights into dynamical interaction mechanisms." Comprehensive review of coupling function estimation across domains. Reviews of Modern Physics 89(4):045001.
• Dörfler, F. & Bullo, F. (2014): "Synchronization in complex networks of phase oscillators: A survey." Derived the critical coupling condition using algebraic connectivity. Automatica 50(6):1539-1564.
• Watts, D. J. & Strogatz, S. H. (1998): "Collective dynamics of 'small-world' networks." The distance-decay model produces topologies resembling small-world networks. Nature 393(6684):440-442.
• Arenas, A. et al. (2008): "Synchronization in complex networks." Comprehensive review of synchronisation on different topologies. Physics Reports 469(3):93-153.

Bayesian Auto-Tune

The prior enables Bayesian optimisation of coupling parameters:

\[p(K_{\text{base}}, \alpha_d | \text{data}) \propto p(\text{data} | K_{\text{base}}, \alpha_d) \cdot p(K_{\text{base}}, \alpha_d)\]

The log_probability provides \(\log p(K_{\text{base}}, \alpha_d)\). When combined with a likelihood (e.g., from observed \(R\) values), this enables posterior inference with only 2 parameters instead of \(N^2\).

Relation to Bayesian Deep Learning

The universal prior plays the same role as weight priors in Bayesian neural networks. Just as a Gaussian prior on neural network weights regularises toward small, general solutions, the coupling prior regularises the SCPN toward coupling topologies that are known to work across domains.

---

3. Pipeline Position

 Domain specification (new domain, no prior knowledge)
                      │
                      ↓
 ┌── UniversalPrior ──────────────────────────────┐
 │                                                 │
 │  .default() → MAP estimate (K=0.47, α=0.25)   │
 │  .sample(rng) → random draw from prior         │
 │  .estimate_Kc(omegas, n) → prior + K_c         │
 │  .log_probability(K, α) → Bayesian score       │
 │                                                 │
 └──────────────────────┬──────────────────────────┘
                        │
                        ↓
        CouplingPrior(K_base, decay_alpha, K_c_estimate)
                        │
                        ↓
        CouplingBuilder.build(n_layers, K_base, decay_alpha)
                        │
                        ↓
        CouplingState(knm, alpha) ──→ UPDEEngine / SplittingEngine

Input Contracts

UniversalPrior constructor:

Parameter	Type	Default	Meaning
K_base_mean	float	0.47	Mean of base coupling prior
K_base_std	float	0.09	Std of base coupling prior
decay_alpha_mean	float	0.25	Mean of distance-decay prior
decay_alpha_std	float	0.07	Std of distance-decay prior

log_probability:

Parameter	Type	Range	Meaning
K_base	float	\(> 0\)	Base coupling to evaluate
decay_alpha	float	\(> 0\)	Distance-decay to evaluate

estimate_Kc:

Parameter	Type	Shape	Meaning
omegas	NDArray[float64]	(N,)	Natural frequencies
n_layers	int	scalar	Number of oscillators/layers

Output Contracts

Method	Returns	Type
default()	MAP coupling config	CouplingPrior
sample(rng)	Random draw	CouplingPrior
estimate_Kc(omegas, n)	Prior + \(K_c\)	CouplingPrior
log_probability(K, α)	Unnormalised log-density	float

---

4. Features

• Empirical prior from 25 cross-domain configurations — not ad hoc, grounded in systematic analysis
• Dimensionality collapse — reduces \(N^2\) coupling parameters to 2
• Bayesian integration — log_probability for posterior inference
• Critical coupling estimate — estimate_Kc combines prior with algebraic connectivity
• Sampling — sample() for Monte Carlo or Bayesian optimisation
• MAP default — default() returns the maximum a posteriori estimate (mean of the prior)
• Rust FFI for log_probability — 2.4x speedup for inner loops of Bayesian optimisation
• Distance-decay matrix — Rust engine builds the \(N \times N\) coupling matrix from \((K_{\text{base}}, \alpha_d)\)
• Domain-agnostic — works for any oscillator network without domain-specific tuning

---

5. Usage Examples

Basic: Get Default Coupling

from scpn_phase_orchestrator.coupling.prior import UniversalPrior

prior = UniversalPrior()
config = prior.default()
print(f"K_base = {config.K_base}")       # 0.47
print(f"decay_alpha = {config.decay_alpha}")  # 0.25

Sample from Prior

import numpy as np
from scpn_phase_orchestrator.coupling.prior import UniversalPrior

prior = UniversalPrior()
rng = np.random.default_rng(42)

for i in range(5):
    config = prior.sample(rng)
    print(f"Sample {i}: K={config.K_base:.3f}, α={config.decay_alpha:.3f}")

For reproducible sampling without passing an explicit generator, use an integer seed in the unsigned 64-bit range:

config = prior.sample(seed=1234)

Boolean, negative, fractional, string, and out-of-range seeds are rejected.

Estimate Critical Coupling

import numpy as np
from scpn_phase_orchestrator.coupling.prior import UniversalPrior

prior = UniversalPrior()
omegas = np.array([1.0, 1.2, 0.8, 1.5, 0.9, 1.1, 1.3, 0.7])
config = prior.estimate_Kc(omegas, n_layers=8)

print(f"K_base = {config.K_base:.3f}")
print(f"K_c = {config.K_c_estimate:.3f}")
print(f"Supercritical: {config.K_base > config.K_c_estimate}")

Bayesian Parameter Sweep

import numpy as np
from scpn_phase_orchestrator.coupling.prior import UniversalPrior

prior = UniversalPrior()

# Grid of (K_base, decay_alpha) values
K_vals = np.linspace(0.1, 1.0, 50)
alpha_vals = np.linspace(0.05, 0.5, 50)

# Log-prior surface
log_prior = np.zeros((50, 50))
for i, K in enumerate(K_vals):
    for j, alpha in enumerate(alpha_vals):
        log_prior[i, j] = prior.log_probability(K, alpha)

# Find MAP
i_max, j_max = np.unravel_index(np.argmax(log_prior), log_prior.shape)
print(f"MAP: K={K_vals[i_max]:.3f}, α={alpha_vals[j_max]:.3f}")

Build Coupling Matrix from Prior

import numpy as np
from scpn_phase_orchestrator.coupling.prior import UniversalPrior
from scpn_phase_orchestrator.coupling.knm import CouplingBuilder

prior = UniversalPrior()
config = prior.default()

cb = CouplingBuilder()
cs = cb.build(
    n_layers=16,
    base_strength=config.K_base,
    decay_alpha=config.decay_alpha,
)
print(f"Coupling matrix shape: {cs.knm.shape}")
print(f"Max coupling: {cs.knm.max():.4f}")
print(f"Min off-diag: {cs.knm[cs.knm > 0].min():.4f}")

---

6. Technical Reference

Class: UniversalPrior

prior

Empirical domain-agnostic prior for coupling hyperparameters.

UniversalPrior provides default/sample/log-probability helpers for K_base and decay_alpha, plus a Dörfler-Bullo-style critical-coupling estimate over the spectral module. When the optional Rust kernel is importable the log-probability path dispatches there; otherwise the NumPy scalar fallback preserves the same Gaussian-prior contract.

Classes

CouplingPrior dataclass

CouplingPrior(
    K_base: float, decay_alpha: float, K_c_estimate: float
)

Coupling configuration: base strength, decay, and K_c estimate.

UniversalPrior

UniversalPrior(
    K_base_mean: float = _K_BASE_MEAN,
    K_base_std: float = _K_BASE_STD,
    decay_alpha_mean: float = _DECAY_ALPHA_MEAN,
    decay_alpha_std: float = _DECAY_ALPHA_STD,
)

Domain-agnostic coupling prior from 25-domainpack empirical distribution.

K_base ~ N(0.47, 0.09), decay_alpha ~ N(0.25, 0.07). Any new domain starts from this prior. Combined with Dörfler-Bullo K_c, collapses auto-tune from 5D optimization to 2D.

Source: R4-A3 cross-domain transfer analysis (Stankovski 2017, Rev. Mod. Phys.).

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def __init__(
    self,
    K_base_mean: float = _K_BASE_MEAN,
    K_base_std: float = _K_BASE_STD,
    decay_alpha_mean: float = _DECAY_ALPHA_MEAN,
    decay_alpha_std: float = _DECAY_ALPHA_STD,
):
    self._K_base_mean = _validate_finite_real(K_base_mean, name="K_base_mean")
    self._K_base_std = _validate_positive_real(K_base_std, name="K_base_std")
    self._decay_alpha_mean = _validate_finite_real(
        decay_alpha_mean, name="decay_alpha_mean"
    )
    self._decay_alpha_std = _validate_positive_real(
        decay_alpha_std, name="decay_alpha_std"
    )

Methods:

sample

sample(
    rng: Generator | None = None, seed: int | None = None
) -> CouplingPrior

Draw a random coupling configuration from the prior.

Pass rng for an explicit generator, or seed to create a seeded one. If neither is given, a fresh unseeded generator is used (NOT reproducible across sessions).

Parameters

rng : np.random.Generator | None NumPy random generator, or None to seed from seed. seed : int | None Seed for the deterministic RNG.

Returns

CouplingPrior A coupling configuration sampled from the prior.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def sample(
    self,
    rng: np.random.Generator | None = None,
    seed: int | None = None,
) -> CouplingPrior:
    """Draw a random coupling configuration from the prior.

    Pass ``rng`` for an explicit generator, or ``seed`` to create a
    seeded one. If neither is given, a fresh unseeded generator is used
    (NOT reproducible across sessions).

    Parameters
    ----------
    rng : np.random.Generator | None
        NumPy random generator, or ``None`` to seed from ``seed``.
    seed : int | None
        Seed for the deterministic RNG.

    Returns
    -------
    CouplingPrior
        A coupling configuration sampled from the prior.
    """
    if rng is None:
        seed = _validate_seed(seed)
        rng = np.random.default_rng(seed)
    K = max(0.01, rng.normal(self._K_base_mean, self._K_base_std))
    alpha = max(0.01, rng.normal(self._decay_alpha_mean, self._decay_alpha_std))
    return CouplingPrior(K_base=K, decay_alpha=alpha, K_c_estimate=0.0)

default

default() -> CouplingPrior

Return the MAP (maximum a posteriori) estimate = the means.

Returns

CouplingPrior Return the MAP (maximum a posteriori) estimate = the means.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def default(self) -> CouplingPrior:
    """Return the MAP (maximum a posteriori) estimate = the means.

    Returns
    -------
    CouplingPrior
        Return the MAP (maximum a posteriori) estimate = the means.
    """
    return CouplingPrior(
        K_base=self._K_base_mean,
        decay_alpha=self._decay_alpha_mean,
        K_c_estimate=0.0,
    )

estimate_Kc

estimate_Kc(
    omegas: FloatArray, n_layers: int
) -> CouplingPrior

Combine prior with Dörfler-Bullo K_c for given omegas.

K_c = max|ω_i - ω_j| / λ₂(L) where L is built from the prior's decay_alpha on a chain graph of n_layers.

Parameters

omegas : FloatArray Natural frequencies in rad/s, shape (N,). n_layers : int Number of SCPN hierarchy layers.

Returns

CouplingPrior The prior combined with the Dörfler-Bullo K_c for the given frequencies.

Raises

TypeError If an argument has the wrong type. ValueError If omegas or the layer count is invalid.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def estimate_Kc(self, omegas: FloatArray, n_layers: int) -> CouplingPrior:
    """Combine prior with Dörfler-Bullo K_c for given omegas.

    K_c = max|ω_i - ω_j| / λ₂(L) where L is built from the prior's
    decay_alpha on a chain graph of n_layers.

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    n_layers : int
        Number of SCPN hierarchy layers.

    Returns
    -------
    CouplingPrior
        The prior combined with the Dörfler-Bullo ``K_c`` for the given frequencies.

    Raises
    ------
    TypeError
        If an argument has the wrong type.
    ValueError
        If ``omegas`` or the layer count is invalid.
    """
    if isinstance(n_layers, bool) or not isinstance(n_layers, Integral):
        raise TypeError(f"n_layers must be an integer, got {n_layers!r}")
    n_layers = int(n_layers)
    if n_layers <= 0:
        raise ValueError("n_layers must be a positive integer")
    from scpn_phase_orchestrator.coupling.spectral import critical_coupling

    prior = self.default()
    omega_values = _validate_frequency_vector(omegas)
    if omega_values.size != n_layers:
        raise ValueError("n_layers must match the length of omegas")
    idx = np.arange(n_layers)
    dist = np.abs(idx[:, np.newaxis] - idx[np.newaxis, :])
    knm: FloatArray = prior.K_base * np.exp(-prior.decay_alpha * dist)
    np.fill_diagonal(knm, 0.0)
    K_c = critical_coupling(omega_values, knm)
    return CouplingPrior(
        K_base=prior.K_base,
        decay_alpha=prior.decay_alpha,
        K_c_estimate=K_c,
    )

log_probability

log_probability(K_base: float, decay_alpha: float) -> float

Log-probability under the Gaussian prior (unnormalised).

Parameters

K_base : float Base coupling strength before spatial modulation. decay_alpha : float Exponential decay rate of the coupling across layer separation.

Returns

float The unnormalised log-probability of the configuration under the prior.

Raises

TypeError If an argument has the wrong type. ValueError If K_base or decay_alpha is out of range.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def log_probability(self, K_base: float, decay_alpha: float) -> float:
    """Log-probability under the Gaussian prior (unnormalised).

    Parameters
    ----------
    K_base : float
        Base coupling strength before spatial modulation.
    decay_alpha : float
        Exponential decay rate of the coupling across layer separation.

    Returns
    -------
    float
        The unnormalised log-probability of the configuration under the prior.

    Raises
    ------
    TypeError
        If an argument has the wrong type.
    ValueError
        If ``K_base`` or ``decay_alpha`` is out of range.
    """
    if isinstance(K_base, bool) or not isinstance(K_base, Real):
        raise TypeError("K_base must be a finite real value")
    if isinstance(decay_alpha, bool) or not isinstance(decay_alpha, Real):
        raise TypeError("decay_alpha must be a finite real value")
    K_base = float(K_base)
    decay_alpha = float(decay_alpha)
    if not np.isfinite(K_base) or not np.isfinite(decay_alpha):
        raise ValueError("K_base and decay_alpha must be finite real values")
    if _HAS_RUST:
        log_prob = float(
            _rust_log_prob(
                K_base,
                decay_alpha,
                self._K_base_mean,
                self._K_base_std,
                self._decay_alpha_mean,
                self._decay_alpha_std,
            )
        )
        if not np.isfinite(log_prob):
            raise ValueError(
                "Rust prior log-probability must be a finite real value"
            )
        return log_prob
    lp_K = -0.5 * ((K_base - self._K_base_mean) / self._K_base_std) ** 2
    lp_a = (
        -0.5 * ((decay_alpha - self._decay_alpha_mean) / self._decay_alpha_std) ** 2
    )
    return lp_K + lp_a

Dataclass: CouplingPrior

@dataclass
class CouplingPrior:
    K_base: float        # Base coupling strength
    decay_alpha: float   # Distance-decay exponent
    K_c_estimate: float  # Critical coupling estimate (0 if not computed)

Empirical Constants

Constant	Value	Meaning	Source
_K_BASE_MEAN	0.47	Mean base coupling	25-domainpack analysis
_K_BASE_STD	0.09	Std base coupling	25-domainpack analysis
_DECAY_ALPHA_MEAN	0.25	Mean distance-decay	25-domainpack analysis
_DECAY_ALPHA_STD	0.07	Std distance-decay	25-domainpack analysis

Rust Engine Functions

// Unnormalised log-probability under bivariate Gaussian
pub fn log_probability(
    k_base: f64, decay_alpha: f64,
    k_mean: f64, k_std: f64,
    alpha_mean: f64, alpha_std: f64,
) -> f64

// Build distance-decay coupling matrix
pub fn distance_decay_matrix(
    n: usize, k_base: f64, decay_alpha: f64,
) -> Vec<f64>   // N×N row-major, diagonal = 0

Auto-Select Logic

try:
    from spo_kernel import prior_log_probability_rust as _rust_log_prob
    _HAS_RUST = True
except ImportError:
    _HAS_RUST = False

Only log_probability uses the Rust path. The sample, default, and estimate_Kc methods are pure Python (no performance benefit from Rust for single calls).

---

7. Performance Benchmarks

Measured on Intel Core i5-11600K @ 3.90 GHz, 32 GB DDR4-2400. Median of 10,000 iterations.

log_probability

Backend	Time (µs)	Speedup
Python	0.60	—
Rust	0.25	2.4x

Why Only 2.4x?

The log_probability computation is 4 floating-point operations (2 subtracts, 2 divides, 2 squares, 1 add). At this scale, the FFI call overhead (~100 ns) is significant relative to the compute time (~150 ns for Rust vs ~600 ns for Python). The speedup is meaningful only in inner loops that call log_probability thousands of times (e.g., MCMC sampling, Bayesian optimisation).

distance_decay_matrix (Rust only)

N	Time (µs)
16	0.8
64	8.5
256	120

This function is Rust-only (no Python fallback exposed via FFI yet). For \(N = 256\), building the coupling matrix takes ~120 µs — negligible compared to integration costs.

Memory Usage

• CouplingPrior dataclass: 3 floats (24 bytes)
• UniversalPrior instance: 4 floats (32 bytes)
• Distance-decay matrix: \(N^2\) floats

Test Coverage

• Rust tests: 7 (prior module in spo-engine)
• Log-probability at mean (= 0), decreases away from mean, symmetric around mean, distance-decay diagonal zero, distance-decay symmetric, distance-decay monotone decreasing, distance-decay values with zero decay
• Python tests: 10 (tests/test_coupling_prior.py)
• Default values, sample reproducibility, sample positive, log_probability at mean, log_probability away, estimate_Kc positive, estimate_Kc with identical omegas, K_c zero for identical omegas, pipeline wiring, Bayesian sweep
• Source lines: 119 (Rust) + 113 (Python) = 232 total

---

8. Citations

1. Stankovski, T., Pereira, T., McClintock, P. V. E., & Stefanovska, A. (2017). "Coupling functions: Universal insights into dynamical interaction mechanisms." Reviews of Modern Physics 89(4):045001. DOI: 10.1103/RevModPhys.89.045001

2. Dörfler, F. & Bullo, F. (2014). "Synchronization in complex networks of phase oscillators: A survey." Automatica 50(6):1539-1564. DOI: 10.1016/j.automatica.2014.04.012

3. Watts, D. J. & Strogatz, S. H. (1998). "Collective dynamics of 'small-world' networks." Nature 393(6684):440-442. DOI: 10.1038/30918

4. Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y., & Zhou, C. (2008). "Synchronization in complex networks." Physics Reports 469(3):93-153. DOI: 10.1016/j.physrep.2008.09.002

5. Fiedler, M. (1973). "Algebraic connectivity of graphs." Czechoslovak Mathematical Journal 23(98):298-305.

6. Barabási, A.-L. & Albert, R. (1999). "Emergence of scaling in random networks." Science 286(5439):509-512. DOI: 10.1126/science.286.5439.509

7. Newman, M. E. J. (2003). "The structure and function of complex networks." SIAM Review 45(2):167-256. DOI: 10.1137/S003614450342480

8. Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence. Springer. ISBN: 978-3-642-69691-6.

---

Edge Cases and Limitations

Very Small \(\sigma_K\) or \(\sigma_\alpha\)

If the standard deviation is set near zero, log_probability returns extreme negative values for any point away from the mean. This effectively makes the prior a delta function. For Bayesian optimisation, \(\sigma \geq 0.01\) is recommended to maintain exploration.

Negative K_base or decay_alpha

The sample() method clamps both parameters to \(\geq 0.01\). Negative base coupling is physically meaningless (repulsive coupling should be modelled differently), and negative decay is unbounded growth with distance.

\(K_c\) Estimation for Disconnected Topologies

If the distance-decay matrix has \(\lambda_2(L) = 0\) (disconnected graph), \(K_c = \infty\). In practice, for \(\alpha_d < 2\) and \(N \geq 4\), the matrix is always connected (all off-diagonal entries are positive, albeit small).

Prior Mismatch

The universal prior is a Gaussian fit to 25 domains. A specific domain may have optimal parameters far from the prior mean. In this case, the prior serves as a starting point, and posterior inference (using log_probability as the prior term) will shift toward the domain-optimal values after a few Bayesian updates.

Independence Assumption

The prior treats \(K_{\text{base}}\) and \(\alpha_d\) as independent. In reality, domains with high base coupling tend to have higher decay (to keep total coupling bounded). A future extension could use a full covariance matrix, but the independent approximation works well for the auto-tune use case.

---

Integration with Other SPO Modules

With CouplingBuilder

The UniversalPrior provides parameters for CouplingBuilder:

prior = UniversalPrior()
config = prior.default()
cs = CouplingBuilder().build(
    n_layers=N,
    base_strength=config.K_base,
    decay_alpha=config.decay_alpha,
)

This is the standard initialisation path for new domains.

With SINDy Auto-Tune

The log_probability serves as a regulariser in SINDy-based coupling estimation:

\[\text{objective} = ||X_{t+1} - f(X_t; K)||^2 + \lambda \cdot (-\log p(K))\]

The prior penalises coupling estimates that deviate from the universal distribution, preventing overfitting to noisy data.

With ActiveInferenceAgent

The prior enables the Active Inference agent to form beliefs about coupling parameters before observing any data. The agent starts with the prior, collects observations, and updates its beliefs using Bayesian inference with log_probability as the prior term.

---

Troubleshooting

Issue: estimate_Kc Returns Very Large K_c

Diagnosis: The algebraic connectivity \(\lambda_2\) is small, meaning the prior's topology is poorly connected. This happens when \(\alpha_d\) is large (steep distance decay).

Solution: Reduce decay_alpha or use a custom prior with smaller \(\alpha_d\) mean.

Issue: sample() Always Returns Similar Values

Diagnosis: The prior is narrow (\(\sigma_K = 0.09\), \(\sigma_\alpha = 0.07\)). 95% of samples fall within \([0.29, 0.65]\) for \(K_{\text{base}}\) and \([0.11, 0.39]\) for \(\alpha_d\).

Solution: For wider exploration, construct a prior with larger standard deviations: UniversalPrior(K_base_std=0.2, decay_alpha_std=0.15).