Universal Coupling Prior

July 2, 2026 · View on GitHub

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:

Scale invariance: The normalised coupling $K/K_c$ has similar optimal values across domains (typically $K/K_c \in [1.5, 3.0]$ ).
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.
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

::: scpn_phase_orchestrator.coupling.prior

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

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
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
Watts, D. J. & Strogatz, S. H. (1998). "Collective dynamics of 'small-world' networks." Nature 393(6684):440-442. DOI: 10.1038/30918
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
Fiedler, M. (1973). "Algebraic connectivity of graphs." Czechoslovak Mathematical Journal 23(98):298-305.
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
Newman, M. E. J. (2003). "The structure and function of complex networks." SIAM Review 45(2):167-256. DOI: 10.1137/S003614450342480
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.