qig.exponential_family

Quantum exponential families and Fisher information geometry.

The core class qig.exponential_family.QuantumExponentialFamily provides several ways to compute the density-matrix derivatives \(\partial\rho / \partial\theta_a\):

method='sld': symmetric logarithmic derivative approximation (fast, ~few-percent error in genuinely quantum, non-commuting cases).
method='duhamel': high-precision Duhamel / Wilcox formula using numerical quadrature over \(s \in [0,1]\) (slower, but serves as a reference).
method='duhamel_spectral' (alias 'duhamel_bch'): uses the spectral/BCH representation of \(H\) to evaluate the Duhamel integral as a matrix function \(f(\mathrm{ad}_H)\) with \(f(z) = (e^z - 1)/z\), avoiding explicit quadrature and matching the Lie-closure discussion in the theory sections.

For small finite-dimensional systems (e.g. the qutrit-pair examples used in the origin paper), the spectral/BCH variant is typically the best choice: it is as accurate as the quadrature-based Duhamel evaluation and more efficient, while remaining faithful to the Kubo-Mori / BKM structure. The 'duhamel' method is retained for validation and for comparison with legacy code paths.

Quantum exponential family and BKM metric interface for the quantum inaccessible game.

This module contains: - local operator bases (Pauli, Gell-Mann, generalised Gell-Mann); - construction of the full operator basis {F_a}; - the QuantumExponentialFamily class providing ρ(θ), log-partition ψ(θ),

Fisher/BKM metric G(θ), and the marginal-entropy constraint.

qig.exponential_family.pauli_basis(site, n_sites)[source]

Create Pauli operator basis for site ‘site’ in an n_sites qubit system.

Returns [sigma_x, sigma_y, sigma_z] tensored with identity on other sites.

Return type:

List[ndarray]

Parameters:

site (int)
n_sites (int)

qig.exponential_family.gell_mann_matrices()[source]

Return the 8 Gell-Mann matrices (generators of SU(3)).

Return type:: List[ndarray]

qig.exponential_family.qutrit_basis(site, n_sites)[source]

Create Gell-Mann operator basis for site ‘site’ in n_sites qutrit system.

Return type:

List[ndarray]

Parameters:

site (int)
n_sites (int)

qig.exponential_family.create_operator_basis(n_sites, d)[source]

Create full operator basis {F_a} for quantum exponential family.

Parameters:

n_sites (int) – Number of sites
d (int) – Local dimension (2 for qubits, 3 for qutrits, d>=2 in general)

Return type:

Tuple[list, list]

class qig.exponential_family.QuantumExponentialFamily(n_sites=None, d=2, n_pairs=None, pair_basis=False)[source]

Bases: object

Quantum exponential family: ρ(θ) = exp(∑ θ_a F_a - ψ(θ))

Parameters:

n_sites (int | None)
d (int)
n_pairs (int | None)
pair_basis (bool)

__init__(n_sites=None, d=2, n_pairs=None, pair_basis=False)[source]

Initialise quantum exponential family.

Parameters:

n_sites (int, optional) – Number of subsystems (for local operator basis)
d (int) – Local dimension (2 for qubits, 3 for qutrits)
n_pairs (int, optional) – Number of entangled pairs (for pair basis)
pair_basis (bool) – If True, use su(d²) operators for entangled pairs If False, use local operators (default)

Notes

Two modes of operation:

Local operators (pair_basis=False, default): - Specify n_sites: number of independent subsystems - Uses operators like σ_x⊗I, I⊗σ_y (local Pauli/Gell-Mann) - Can ONLY represent separable states - Suitable for testing, but NOT for quantum game with entanglement
Pair operators (pair_basis=True): - Specify n_pairs: number of entangled pairs - Uses su(d²) generators acting on each pair - Can represent entangled states (including Bell states) - Required for proper quantum inaccessible game - Fisher metric G is block-diagonal (one block per pair)

rho_from_theta(theta)[source]

Compute ρ(θ) = exp(K(θ) - ψ(θ)) where K(θ) = ∑ θ_a F_a.

Return type:: ndarray
Parameters:: theta (ndarray)

psi(theta)[source]

Compute the cumulant generating function ψ(θ) = log Tr(exp(∑ θ_a F_a)).

This is the log partition function for the exponential family.

Return type:: float
Parameters:: theta (ndarray)

log_partition(theta)[source]

Deprecated: Use psi() instead.

Compute the cumulant generating function ψ(θ) = log Tr(exp(∑ θ_a F_a)).

Return type:: float
Parameters:: theta (ndarray)

rho_derivative(theta, a, **kwargs)[source]

Compute ∂ρ/∂θ_a using the correct quantum formula.

Three main methods available. 1. ‘sld’: Symmetric Logarithmic Derivative (fast, ~5% error)

∂ρ/∂θ = (1/2)[ρ(F - ⟨F⟩I) + (F - ⟨F⟩I)ρ]

‘duhamel’: Duhamel exponential formula with numerical quadrature (slower, ~5e-6 error for n_points≈200)

∂ρ/∂θ = ∫₀¹ exp(sH)(F - ⟨F⟩I)exp((1-s)H) ds
‘duhamel_spectral’ / ‘duhamel_bch’: Duhamel formula via spectral/BCH representation of H (no explicit s-quadrature, high precision).
‘duhamel_block’: Higham block-matrix identity for the Fréchet derivative (one 2D×2D matrix exponential; no eigendecomposition, no quadrature).

⚠️ QUANTUM ALERT: Simple ρ(F - ⟨F⟩I) is WRONG for non-commuting operators!

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
a (int) – Parameter index
method (str, default 'sld') – ‘sld’ for fast SLD, ‘duhamel’ for high precision
n_points (int, default 100) – Quadrature points for Duhamel (ignored for ‘sld’)

Returns:

drho – Derivative ∂ρ/∂θ_a (Hermitian matrix)

Return type:

ndarray, shape (D, D)

Notes

Quantum derivative principles applied: ✅ Check operator commutation: Uses symmetric/integral form ✅ Verify operator ordering: Preserves ordering in exponentials ✅ Distinguish quantum vs classical: Uses quantum formulas ✅ Respect Hilbert space structure: Preserves Hermiticity

rho_second_derivative(theta, a, b, method='numerical_duhamel', n_points=100, eps=1e-07)[source]

Compute ∂²ρ/∂θ_a∂θ_b using high-precision numerical differentiation.

Strategy: Use finite differences of the high-precision Duhamel ∂ρ/∂θ.

∂²ρ/∂θ_a∂θ_b ≈ [∂ρ/∂θ_a(θ+ε·e_b) - ∂ρ/∂θ_a(θ-ε·e_b)] / (2ε)

This gives 0.55-2.6% error, which is 30-100× better than the analytic SLD-based formula.

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
a (int) – First parameter index
b (int) – Second parameter index
method (str, default 'numerical_duhamel') – Currently only ‘numerical_duhamel’ is supported
n_points (int, default 100) – Quadrature points for Duhamel integration
eps (float, default 1e-7) – Finite difference step size

Returns:

d2rho – Second derivative ∂²ρ/∂θ_a∂θ_b (Hermitian matrix)

Return type:

ndarray, shape (D, D)

Notes

Uses central differences for better accuracy
Hermiticity preserved to machine precision
More accurate than analytic SLD-based second derivative

Quantum derivative principles applied: ✅ Respects non-commutativity through Duhamel integration ✅ Preserves Hermiticity ✅ Uses high-precision first derivatives

fisher_information(theta)[source]

Compute Fisher information (BKM metric) G(θ) = ∇∇ψ(θ) using the Kubo-Mori / BKM inner product.

For a quantum exponential family ρ(θ) = exp(K(θ)) / Z(θ) with K(θ) = ∑_a θ_a F_a, the Bogoliubov-Kubo-Mori metric is:

\[G_{ab}(θ) = \int_0^1 \mathrm{Tr}\left( ρ(θ)^s \tilde{F}_a ρ(θ)^{1-s} \tilde{F}_b \right) \mathrm{d}s\]

where F̃_a = F_a - Tr[ρ(θ) F_a] I are centred sufficient statistics.

In the eigenbasis of ρ(θ), this reduces to the spectral representation with the Morozova-Chentsov function c(λ, μ) = (log λ - log μ)/(λ - μ) (diagonal limit: c(λ, λ) = 1/λ).

When all F_a commute with ρ(θ) (classical case), this reduces to the usual covariance Fisher information matrix.

This implementation:

Diagonalises ρ(θ) = U diag(p) U†
Centres each F_a in that basis
Applies the BKM kernel c(p_i, p_j) to all matrix elements
Symmetrises G to guard against numerical asymmetries

Return type:: ndarray
Parameters:: theta (ndarray)

fisher_information_product(theta, check_product=True, tol=1e-06)[source]

Compute Fisher information exploiting block-diagonal structure for product states.

For product states ρ = ρ₁ ⊗ ρ₂ ⊗ … ⊗ ρₙ with pair-based operators, the BKM Fisher metric is block-diagonal:

G = diag(G₁, G₂, …, Gₙ)

where each Gₖ is computed from the marginal ρₖ on pair k.

Complexity:

Full: O(D³) = O(d^(6n)) where D = d^(2n)
Block: O(n × d¹²) - exponentially faster for n > 1

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
check_product (bool, default True) – If True, verify state is approximately a product state. If False, assume product structure (faster but may be wrong).
tol (float, default 1e-6) – Tolerance for product state check

Returns:

G – Block-diagonal Fisher information matrix

Return type:

ndarray, shape (n_params, n_params)

Raises:

ValueError – If not using pair_basis mode If check_product=True and state is not a product state

Notes

For n qutrit pairs (d=3): - n=2: Full O(81³)=530k, Block O(2×3¹²)=1M → similar - n=3: Full O(729³)=387M, Block O(3×3¹²)=1.6M → 240× faster - n=4: Full O(6561³)=282B, Block O(4×3¹²)=2.1M → 134000× faster

The crossover point where block computation wins is typically n≥2 for qutrits and n≥3 for qubits.

marginal_entropy_constraint_theta_only(theta)[source]

Compute constraint C(θ) = ∑ᵢ hᵢ and gradient ∇C using θ-only formulas.

This is a fast, exact method that avoids materializing ∂ρ/∂θ for each parameter. Instead, it uses the BKM inner product:

∂C/∂θ_a = ⟨F̃_a, B⟩_BKM

where:: F̃_a = F_a - ⟨F_a⟩I (centered operator) B = ∑ᵢ Bᵢ (lifted test operator) Bᵢ = (log ρᵢ + Iᵢ) ⊗ I_rest

This reuses the eigendecomposition and BKM kernel from fisher_information(), achieving machine precision (~10⁻¹⁴) with ~100× speedup over Duhamel method.

Parameters:

theta (ndarray) – Natural parameters

Return type:

Tuple[float, ndarray]

Returns:

C (float) – Constraint value ∑ᵢ hᵢ
grad_C (ndarray, shape (n_params,)) – Gradient ∇C

marginal_entropy_constraint(theta, method='theta_only')[source]

Compute constraint value C(θ) = ∑_i h_i and gradient ∇C.

This method dispatches to different implementations based on the method parameter.

Parameters:

theta (ndarray) – Natural parameters
method (str, default 'theta_only') –
Method for gradient computation. Options:
- 'theta_only': Fast θ-only method using BKM inner products (default). ~100× faster, machine precision accuracy.
- 'duhamel': Legacy method materializing ∂ρ/∂θ via Duhamel. Slow but kept for verification.
- 'sld': Legacy method using SLD approximation. Faster than Duhamel but ~5% error.

Return type:

Tuple[float, ndarray]

Returns:

C (float) – Constraint value ∑ᵢ hᵢ
grad_C (ndarray) – Gradient ∇C

third_cumulant_contraction(theta, method='fd')[source]

Compute (∇G)[θ], the third cumulant tensor contracted with θ.

This is the matrix with (i,j) entry: ∑_k (∂G_ik/∂θ_j) θ_k

Following the paper’s Appendix (eq. 824-826), this appears in the Jacobian as:: M = -G - (∇G)[θ] + …

The third cumulant ∇G = ∇³ψ is totally symmetric in all three indices.

We compute ∂G_ab/∂θ_c by differentiating the spectral BKM formula using perturbation theory for eigenvalues and eigenvectors.

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
method (str, optional) – Method for computing the third cumulant: - ‘fd’ (default): Finite differences of Fisher metric (fast, accurate) - ‘analytic’: Analytic perturbation theory (slow, approximate)

Returns:

contraction – Matrix (∇G)[θ] with (i,j) entry = ∑_k (∂G_ik/∂θ_j) θ_k

Return type:

ndarray, shape (n_params, n_params)

Notes

Quantum derivative principles applied: ✅ Check operator commutation: F_a, F_b, F_c may not commute ✅ Verify operator ordering: Careful in spectral differentiation ✅ Distinguish quantum vs classical: Uses quantum covariance derivatives ✅ Respect Hilbert space structure: Works on full Hilbert space ✅ Question each derivative step: Uses perturbation theory

The finite difference method is much faster (~100-500×) and avoids expensive ∂ρ/∂θ computations.

psi_hessian_block(theta, param_indices=None)[source]

Compute Hessian of ψ(θ) = log Tr(exp(K(θ))) via block Fréchet derivatives.

Return type:

ndarray

Parameters:

theta (ndarray)
param_indices (List[int] | None)

This returns the (selected) 2nd cumulant matrix:: ∂²ψ/∂θ_a∂θ_b

Notes

This method is intended for validation / small systems. It scales as O(m²) block matrix exponentials, where m = len(param_indices) (or n_params).

constraint_hessian_fd_theta_only(theta, eps=1e-05)[source]

Compute constraint Hessian ∇²C using finite differences of θ-only gradient.

This is a fast, accurate method that computes:: ∂²C/∂θ_a∂θ_b ≈ [∇C(θ + eps·e_b) - ∇C(θ - eps·e_b)]_a / (2·eps)

By differentiating the exact θ-only gradient, this achieves better accuracy than the current approach which uses exact formulas with approximate second derivatives (FD of Duhamel drho).

Expected speedup: 50-100× over current Duhamel-based method. Expected accuracy: ~10⁻⁸ (better than current ~10⁻⁶).

Parameters:

theta (ndarray) – Natural parameters
eps (float, optional) – Finite difference step size (default: 1e-5) Optimal for central differences: h ≈ (machine_eps)^(1/3) ≈ 1e-5

Returns:

hess – Hessian matrix ∇²C, symmetric real matrix

Return type:

ndarray, shape (n_params, n_params)

Notes

Why this is better than current approach:

Current (two approximations): Error ≈ 10⁻⁶

∂²C/∂θ_a∂θ_b = f(∂²ρ/∂θ_a∂θ_b) = f(FD(∂ρ/∂θ)) = f(FD(Duhamel(ρ)))

New (one approximation): Error ≈ 10⁻⁸

∂²C/∂θ_a∂θ_b ≈ FD(∂C/∂θ_a) = FD(exact_BKM_formula(ρ))

Key insight: Differentiating an exact gradient is more accurate than using an exact formula with approximate second derivatives.

constraint_hessian(theta, method='fd_theta_only', n_points=100, eps=1e-05)[source]

Compute ∇²C, the Hessian of the constraint C(θ) = ∑ᵢ hᵢ(θ).

This method dispatches to different implementations based on the method parameter.

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
method (str, default 'fd_theta_only') –
Method for Hessian computation. Options:
- 'fd_theta_only': FD of θ-only gradient (default). ~50-100× faster, ~10⁻⁸ error, recommended.
- 'duhamel': FD of Duhamel drho (legacy, slow, ~10⁻⁶ error).
- 'sld': Analytic formula using SLD (legacy, fast but ~10% error).
n_points (int, default 100) – Quadrature points for Duhamel (ignored for other methods)
eps (float, default 1e-5) – Finite difference step size

Returns:

hessian – Constraint Hessian ∇²C (symmetric matrix)

Return type:

ndarray, shape (n_params, n_params)

Notes

The default ‘fd_theta_only’ method uses finite differences of the exact θ-only gradient, achieving better accuracy and much better performance than the legacy methods which use exact formulas with approximate inputs.

lagrange_multiplier_gradient(theta, method='sld', n_points=100)[source]

Compute ∇ν, the gradient of the Lagrange multiplier ν(θ).

From the paper (equations 831-832), the gradient components are:

∂ν/∂θ_j = (1/||a||²) [a^T G e_j + a^T (∇G)[θ] e_j + (∇a)_j^T G θ - 2ν a^T (∇a)_j]

where:

ν = (a^T G θ)/(a^T a) is the Lagrange multiplier
a = ∇C = ∇(∑h_i) is the constraint gradient
G = Fisher information (BKM metric)
(∇G)[θ] = third cumulant tensor contracted with θ
∇a = ∇²C is the constraint Hessian

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
method (str, default 'sld') – ‘sld’ or ‘duhamel’ for derivative precision
n_points (int, default 100) – Quadrature points for Duhamel (ignored for ‘sld’)

Returns:

grad_nu – Gradient ∇ν

Return type:

ndarray, shape (n_params,)

Notes

Uses all the high-precision components from Steps 1-3: - BKM metric G (spectral formula) - Third cumulant (∇G)[θ] (perturbation theory) - Constraint Hessian ∇²C (Duhamel for high precision)

jacobian(theta, method='duhamel', n_points=100)[source]

Compute the Jacobian M = ∂F/∂θ of the constrained dynamics.

From the paper (equations 824-827):

F(θ) = -G(θ)θ + ν(θ)a(θ) M = ∂F/∂θ = -G - (∇G)[θ] + ν∇²C + a(∇ν)^T

For systems with local operators only (no entanglement):

Structural identity Gθ = -a holds
ν = -1 always, ∇ν = 0
Simplifies to: M = -G - (∇G)[θ] - ∇²C

For systems with entangling operators (pair basis):

Structural identity BROKEN: Gθ ≠ -a
ν ≠ -1, ∇ν ≠ 0
Must use full formula: M = -G - (∇G)[θ] + ν∇²C + a(∇ν)^T

Parameters:

theta (ndarray, shape (n_params,)) – Natural parameters
method (str, default 'duhamel') – ‘sld’ or ‘duhamel’ for derivative precision Changed default to ‘duhamel’ for better accuracy with entangled states
n_points (int, default 100) – Quadrature points for Duhamel (ignored for ‘sld’)

Returns:

M – Jacobian matrix

Return type:

ndarray, shape (n_params, n_params)

Notes

This is the full Jacobian for GENERIC dynamics:: θ̇ = F(θ) = -Gθ + νa

The degeneracy of M determines the geometry of the constraint manifold and is central to the paper’s analysis of the inaccessible game.

symmetric_part(theta, method='duhamel', n_points=100)[source]

Compute symmetric part S of the flow Jacobian M.

The GENERIC decomposition splits M into symmetric and antisymmetric parts: M = S + A, where S = (M + M^T)/2

The symmetric part S generates the irreversible (dissipative) dynamics.

Parameters:

theta (np.ndarray) – Natural parameters
method (str) – Method for computing Jacobian
n_points (int) – Number of points for numerical integration

Returns:

S – Symmetric part of Jacobian, satisfies S = S^T

Return type:

np.ndarray, shape (n_params, n_params)

Notes

The symmetric part satisfies key degeneracy conditions: - S @ a ≈ 0, where a = ∇C is the constraint gradient - θ^T S θ ≥ 0 (entropy production non-negative on tangent space)

See also

symmetric_part: Symmetric part of Jacobian
verify_degeneracy_conditions: Verify GENERIC structure

verify_degeneracy_conditions(theta, method='duhamel', n_points=100, tol=1e-06)[source]

Verify degeneracy conditions for GENERIC structure.

The GENERIC structure requires that the symmetric and antisymmetric parts satisfy specific degeneracy conditions related to the constraint and entropy gradients.

Parameters:

theta (np.ndarray) – Natural parameters
method (str) – Method for computing Jacobian
n_points (int) – Number of points for numerical integration
tol (float) – Tolerance for degeneracy conditions

Returns:

diagnostics – Dictionary containing: - ‘S’: Symmetric part - ‘A’: Antisymmetric part - ‘constraint_gradient’: a = ∇C - ‘entropy_gradient’: ∇H = -G @ theta - ‘S_annihilates_constraint’: ||S @ a|| - ‘A_annihilates_entropy_gradient’: ||A @ (-G @ theta)|| - ‘entropy_production’: θ^T S θ - ‘S_symmetric_error’: ||S - S^T|| - ‘A_antisymmetric_error’: ||A + A^T|| - ‘reconstruction_error’: ||M - (S + A)|| - ‘all_passed’: bool indicating if all checks passed

Return type:

dict

Notes

Degeneracy conditions (should hold within tolerance): 1. S @ a ≈ 0: Symmetric part annihilates constraint gradient 2. A @ ∇H ≈ 0: Antisymmetric part annihilates entropy gradient 3. θ^T S θ ≥ 0: Entropy production non-negative 4. S = S^T: Symmetry (machine precision) 5. A = -A^T: Antisymmetry (machine precision) 6. M = S + A: Reconstruction (machine precision)

von_neumann_entropy(theta)[source]

Compute the von Neumann entropy.

For exponential families ρ(θ) = exp(K(θ) - ψ(θ)), we use the fundamental identity: H(ρ) = ψ(θ) - θ^T ∇ψ(θ)

where ∇ψ(θ) is the analytical gradient of the cumulant generating function ψ(θ).

This directly uses the exponential family structure and is more numerically stable than eigendecomposition for states within the exponential family manifold.

Parameters:: theta (ndarray) – Natural parameters
Returns:: Von Neumann entropy in nats
Return type:: float

mutual_information(theta)[source]

Compute mutual information I = C - H where C = ∑h_i, H = S(ρ).

For separable states: I = 0 For entangled states: I > 0 Maximum for Bell states: I = 2log(d)

Parameters:: theta (ndarray) – Natural parameters
Returns:: Mutual information in nats
Return type:: float

Notes

This only makes sense for pair-based systems. For local operators, I ≈ 0 always since those can only create separable states.

purity(theta)[source]

Compute purity Tr(ρ²).

Pure states: Tr(ρ²) = 1 Maximally mixed: Tr(ρ²) = 1/D

Parameters:: theta (ndarray) – Natural parameters
Returns:: Purity, between 1/D and 1
Return type:: float

validate_sigma(sigma)[source]

Validate that σ is a valid density matrix of correct dimension.

Parameters:

sigma (ndarray) – Matrix to validate

Return type:

Tuple[bool, str]

Returns:

is_valid (bool) – True if σ is a valid density matrix
message (str) – “Valid” or error description

is_product_sigma(sigma, tol=1e-06)[source]

Check if σ has product structure σ = σ₁⊗σ₂⊗…⊗σₙ.

Parameters:

sigma (ndarray) – Density matrix to check
tol (float) – Tolerance for product structure detection

Return type:

Tuple[bool, Optional[List[ndarray]]]

Returns:

is_product (bool) – True if σ is (approximately) a product state
factors (list of ndarray or None) – If product, the individual σₖ matrices; otherwise None

detect_sigma_structure(sigma)[source]

Detect structure of σ for efficiency optimisation.

Parameters:: sigma (ndarray) – Regularisation matrix
Returns:: structure – One of: ‘isotropic’, ‘product’, ‘pure’, ‘general’
Return type:: str

regularise_pure_state(psi, epsilon, sigma=None)[source]

Create regularised density matrix from pure state.

ρ_ε = (1-ε)|ψ⟩⟨ψ| + ε σ

Parameters:

psi (ndarray) – Pure state vector (length D)
epsilon (float) – Regularisation strength, must be in (0, 1)
sigma (ndarray, optional) – Regularisation direction. Default: I/D (isotropic)

Returns:

rho_epsilon – Regularised density matrix (D×D)

Return type:

ndarray

get_bell_state_parameters(epsilon=1e-06, log_epsilon=None, sigma=None, sigma_per_pair=None, bell_indices=None)[source]

Get natural parameters θ corresponding to a regularised Bell state.

A Bell state is a pure state (rank 1), which lies at the boundary of the exponential family where natural parameters θ → -∞.

For regularised state: ρ_ε = (1-ε)|Φ⟩⟨Φ| + ε σ

We compute θ by solving: ρ_ε = exp(Σ θₐFₐ - ψ(θ))

This gives: θₐ ∝ Tr(log(ρ_ε) Fₐ)

Parameters:

epsilon (float, default 1e-6) – Regularisation parameter. Smaller epsilon → closer to pure Bell state (more negative θ). Must be > 0 to avoid singularities.
log_epsilon (float, optional) – If provided, overrides epsilon via log ε = log_epsilon. This is numerically convenient when exploring very small ε.
sigma (ndarray, optional) – Full D×D regularisation matrix. Any valid density matrix. Use for entangled σ (inter-pair correlations in perturbation). If None and sigma_per_pair is None, uses I/D (isotropic).
bell_indices (list of int, optional) – Which Bell state (k=0,…,d-1) to use for each pair. Default: all zeros (standard Bell state |Φ₀⟩). Allows exploring different pure state origins.
sigma_per_pair (list of ndarray, optional) – List of n density matrices, each d²×d² for one pair. Constructs σ = σ₁⊗σ₂⊗…⊗σₙ (product structure). Efficient O(n) computation preserved. Only valid for multi-pair systems (n_pairs > 1).

Returns:

theta – Natural parameters for the regularized Bell state. Many components will be large and negative (approaching -∞ for pure state).

Return type:

ndarray

Raises:

ValueError – If both sigma and sigma_per_pair are provided. If sigma_per_pair has wrong length or invalid matrices.

Notes

The regularisation matrix σ encodes the “direction of approach” to the pure-state boundary (see CIP-0008 and entropy_time_paths.ipynb):

Different σ = different “meridians” from the north pole
Isotropic σ = I/D gives the “boring” symmetric departure
Anisotropic σ reveals the tangent cone of possible departures

For the exponential family ρ(θ) = exp(H(θ))/Z where H = Σ θₐFₐ, we have: log(ρ) = H - log(Z)·I

Since Tr(Fₐ) = 0 for our operators, we get: θₐ = Tr(log(ρ) Fₐ) / Tr(FₐFₐ)

Multi-pair note: For n_pairs > 1, our operator basis is the direct sum of per-pair su(d²) algebras, NOT the full su(D) algebra. The regularised Bell state may have components outside this subspace. The returned θ is the projection onto our subspace, which is correct for the inaccessible game dynamics (which stay in this subspace). Reconstruction via rho_from_theta(θ) will NOT exactly match the target ρ_ε, but the dynamics are correct.