GENERIC Decomposition

The GENERIC (General Equation for Non-Equilibrium Reversible-Irreversible Coupling) decomposition separates quantum dynamics into reversible (Hamiltonian) and irreversible (dissipative) components.

For the quantum inaccessible game, this decomposition reveals the structure of constrained entropy-maximizing dynamics on the marginal-entropy manifold.

Overview

The flow Jacobian \(M = \partial F/\partial\theta\) naturally decomposes as:

\[M = S + A\]

where:

\(S = \frac{1}{2}(M + M^T)\) is the symmetric part (dissipative/irreversible)
\(A = \frac{1}{2}(M - M^T)\) is the antisymmetric part (reversible/Hamiltonian)

From this decomposition, we can extract:

Effective Hamiltonian \(H_{\text{eff}}\) from \(A\)
Diffusion operator \(\mathcal{D}[\rho]\) from \(S\)

The full dynamics then take the form:

\[\dot{\rho} = -i[H_{\text{eff}}, \rho] + \mathcal{D}[\rho]\]

Quick Start

Basic Example

import numpy as np
from qig.exponential_family import QuantumExponentialFamily
from qig.generic_decomposition import run_generic_decomposition

# Initialize 2-qubit entangled pair
exp_fam = QuantumExponentialFamily(n_pairs=1, d=2, pair_basis=True)

# Choose a state (near LME origin)
theta = 0.1 * np.random.randn(exp_fam.n_params)

# Run complete GENERIC decomposition
results = run_generic_decomposition(
    theta, exp_fam,
    verbose=True,
    print_summary=True
)

This executes all 12 steps and provides comprehensive diagnostics.

Understanding the Results

The results dictionary contains:

# Information geometry
results['psi']      # Cumulant generating function ψ(θ)
results['mu']       # Mean parameters ∇ψ
results['G']        # Fisher information (BKM metric)
results['a']        # Constraint gradient
results['nu']       # Lagrange multiplier

# GENERIC decomposition
results['M']        # Full Jacobian
results['S']        # Symmetric part (dissipative)
results['A']        # Antisymmetric part (reversible)

# Physical operators
results['H_eff']    # Effective Hamiltonian
results['eta']      # Hamiltonian coefficients
results['D_rho']    # Diffusion operator (if computed)

# Diagnostics
results['diagnostics']  # Validation checks

Step-by-Step Guide

Step 1: Initialize Exponential Family

Choose your system:

# For qubits (d=2)
exp_fam = QuantumExponentialFamily(n_pairs=1, d=2, pair_basis=True)

# For qutrits (d=3)
exp_fam = QuantumExponentialFamily(n_pairs=1, d=3, pair_basis=True)

# For multiple pairs
exp_fam = QuantumExponentialFamily(n_pairs=2, d=2, pair_basis=True)

Step 2: Choose Initial State

Select natural parameters \(\theta\):

# At LME origin (maximally entangled)
theta = np.zeros(exp_fam.n_params)

# Near origin (slightly entangled)
theta = 0.1 * np.random.randn(exp_fam.n_params)

# Random state
theta = np.random.randn(exp_fam.n_params)

Step 3: Run Decomposition

Use the convenience function:

results = run_generic_decomposition(
    theta, exp_fam,
    method='duhamel',          # or 'sld' for faster approximation
    compute_diffusion=False,   # Set True to compute D[ρ] (expensive!)
    verbose=True,              # Print progress
    print_summary=True         # Show summary at end
)

Or use the class for more control:

from qig.generic_decomposition import GenericDecomposition

decomp = GenericDecomposition(exp_fam, compute_diffusion=False)
results = decomp.compute_all(theta, verbose=True)
decomp.print_summary(detailed=True)

Step 4: Extract Key Components

Access the decomposition:

# Jacobian decomposition
M = results['M']  # Full Jacobian
S = results['S']  # Symmetric (dissipative)
A = results['A']  # Antisymmetric (reversible)

# Verify: M = S + A
assert np.allclose(M, S + A)

# Effective Hamiltonian
H_eff = results['H_eff']
eta = results['eta']  # Coefficients

# Check Hermiticity
assert np.allclose(H_eff, H_eff.conj().T)

Step 5: Validate Results

Check diagnostics:

diag = results['diagnostics']

# Overall pass/fail
if diag['all_checks_pass']:
    print("✓ All validation checks passed!")
else:
    print("✗ Some checks failed:")
    for name, passed in diag['checks'].items():
        if not passed:
            print(f"  - {name}")

# Detailed error metrics
print(f"S symmetry error: {diag['S_symmetry_error']:.2e}")
print(f"H_eff Hermiticity: {diag['H_eff_hermiticity_error']:.2e}")

Advanced Usage

Computing the Diffusion Operator

The diffusion operator \(\mathcal{D}[\rho]\) maps the symmetric flow to density matrix space. Warning: This is computationally expensive as it requires Kubo-Mori derivatives.

results = run_generic_decomposition(
    theta, exp_fam,
    compute_diffusion=True,  # Enable D[ρ] computation
    verbose=False
)

D_rho = results['D_rho']

# Properties of D[ρ]
assert np.allclose(D_rho, D_rho.conj().T)  # Hermitian
assert abs(np.trace(D_rho)) < 1e-10         # Traceless

Integrating GENERIC Dynamics

Use qig.dynamics.GenericDynamics to track GENERIC structure along trajectories:

from qig.dynamics import GenericDynamics

dyn = GenericDynamics(exp_fam)

# Integrate with monitoring
result = dyn.integrate_with_monitoring(
    theta_0, (0.0, 1.0), n_points=50,
    compute_diffusion=False
)

# Access GENERIC structure along trajectory
H_eff_traj = result['H_eff']           # Hamiltonians
entropy_prod = result['entropy_production']  # dS/dt
S_norms = result['S_norm']             # ||S|| over time
A_norms = result['A_norm']             # ||A|| over time

Separating Reversible and Irreversible Dynamics

Integrate only one component:

# Reversible (Hamiltonian) only
rev_result = dyn.integrate_reversible(theta_0, (0.0, 1.0))

# Irreversible (dissipative) only
irr_result = dyn.integrate_irreversible(theta_0, (0.0, 1.0))

# Compare with full dynamics
full_result = dyn.integrate(theta_0, (0.0, 1.0))

Performance Considerations

Computational Cost

The main computational bottlenecks are:

Fisher information \(G(\theta)\) - \(O(n^2 D^2)\)
Jacobian \(M\) - \(O(n^3)\) for third cumulants
Diffusion operator \(\mathcal{D}[\rho]\) - \(O(n D^4)\) (very expensive!)

For large systems:

Skip diffusion operator computation (compute_diffusion=False)
Use method='sld' for faster (but less accurate) derivatives
Cache structure constants (computed once per basis)

Accuracy vs Speed

Two methods for derivatives:

# High accuracy (~10⁻⁶ error), slow
results = run_generic_decomposition(theta, exp_fam, method='duhamel')

# Fast (~5% error), much faster
results = run_generic_decomposition(theta, exp_fam, method='sld')

Choose based on your needs:

Use 'duhamel' for publication-quality results
Use 'sld' for exploration and prototyping

Troubleshooting

Common Issues

Degeneracy conditions fail

The conditions \(S \cdot a \approx 0\) and \(A \cdot (-G\theta) \approx 0\) may not hold exactly far from equilibrium. This is expected and doesn’t indicate an error in the computation.

Large constraint gradient norm

If \(\|a\|\) is very large, you may be far from the constraint manifold. Project back onto the constraint or reduce parameter magnitudes.

Numerical instabilities

Near pure states or boundaries, numerical precision may degrade. Use:

Smaller parameter values
Higher precision (method='duhamel')
Regularization in structure constant computation

Interpreting Diagnostics

Always passing (algebraic properties):

S_symmetric: \(S = S^T\)
A_antisymmetric: \(A = -A^T\)
M_reconstructs: \(M = S + A\)
H_eff_hermitian: \(H_{\text{eff}} = H_{\text{eff}}^\dagger\)
H_eff_traceless: \(\text{Tr}(H_{\text{eff}}) = 0\)

May fail (physical conditions):

degeneracy_S: \(\|S \cdot a\| < \epsilon\) (depends on state)
degeneracy_A: \(\|A \cdot (-G\theta)\| < \epsilon\) (depends on state)
tangency: Flow tangent to constraint (numerical precision)

Examples

See the examples/ directory for complete examples:

generic_decomposition_demo.py - Basic decomposition
generic_decomposition_complete.py - Full workflow with visualizations