Symbolic Computation for Qutrits

This document describes the symbolic computation approach for deriving analytic expressions in the quantum inaccessible game for qutrit systems.

Overview

The qig.symbolic module provides exact symbolic expressions for the GENERIC decomposition of qutrit pair systems. A key breakthrough is that for maximally entangled states, we can compute exp(K) exactly with no Taylor approximation.

Key results:

  • The antisymmetric part A ≠ 0 for entangled states

  • Exact exp(K) via block decomposition: 9×9 → 3×3 + 2×2 + 1×1×4

  • 20 block-preserving generators span the full entangled subspace

  • All eigenvalues are at most quadratic (no cubic equations)

Natural Parameter Interpretation

The natural parameters θ in the quantum exponential family have specific meanings:

  • θ = 0: Maximally mixed state (ρ = I/D)

  • θ → -∞ (large negative): Locally maximally entangled (LME/Bell) states

For LME states, regularization (ε ~ 10⁻³) keeps parameters finite but large: ||θ|| ~ 3 with some components around -3.

The su(9) Pair Basis

For a pair of qutrits (d=3), we use the full su(9) Lie algebra:

  • 80 generators* (compared to 16 for local su(3)⊗su(3) basis)

  • Can represent entangled states including Bell states

  • Structural identity does not hold: \(G\theta \neq -a\) (unlike separable states)

This breaking of the structural identity is what allows \(A \neq 0\).

LME Block Decomposition (Key Breakthrough)

For locally maximally entangled (LME) states, the full 9×9 eigenvalue problem decomposes into smaller blocks, all with at most quadratic eigenvalues.

The Entangled Subspace

LME states like \(|\psi\rangle = \frac{1}{\sqrt{3}}(|00\rangle + |11\rangle + |22\rangle)\) live in the 3-dimensional subspace \(\{|00\rangle, |11\rangle, |22\rangle\}\).

In a reordered basis, the 9×9 K matrix becomes block diagonal:

\[\begin{split}K_{\text{block}} = \begin{pmatrix} K_{3\times 3} & 0 \\ 0 & K_{6\times 6} \end{pmatrix}\end{split}\]

The 6×6 block further decomposes into 2×2 + 1×1×4.

Block-Preserving Generators

20 generators preserve this block structure:

  • 4 local: \(\lambda_3 \otimes I\), \(I \otimes \lambda_3\), \(\lambda_8 \otimes I\), \(I \otimes \lambda_8\)

  • 16 entangling: \(\lambda_1 \otimes \lambda_1\), \(\lambda_1 \otimes \lambda_2\), etc.

These span the full entangled subspace (rank 9), enabling exploration of all maximally entangled states while maintaining exact computation.

Reduced Density Matrix Structure

Reduced density matrices also have a special block form:

\[\begin{split}\rho_1 = \begin{pmatrix} a & b & 0 \\ b^* & c & 0 \\ 0 & 0 & d \end{pmatrix}\end{split}\]

Lie-Algebraic Origin

This structure comes from SU(3)’s canonical subgroup decomposition.

\[\mathfrak{su}(3) = \mathfrak{su}(2) \oplus \mathfrak{u}(1) \oplus \text{(ladder operators)}\]

The 8 Gell-Mann matrices \(\lambda_1, \ldots, \lambda_8\) split as:

Type

Generators

Action

SU(2)

\(\lambda_1,2,3\)

Mix \(|0\rangle, |1\rangle\)

U(1)

\(\lambda_8\)

Diagonal; separates (0,1) from 2

Ladder

\(\lambda_4,5,6,7\)

Mix \(|2\rangle\) with others

Key insight: Any state lacking coherences with \(|2\rangle\) has zero coefficients for the ladder operators, automatically giving the 2×2 + 1 block form.

For the quantum inaccessible game with Gell-Mann basis, this block structure applies to the reduced density matrices obtained by partial trace. Specifically:

  • The maximally mixed state: \(\rho = I/3\) (corresponds to θ = 0)

  • LME (Bell) states: \(|\Phi^+\rangle\) (corresponds to θ → -∞)

  • Partial traces of maximally entangled states

  • States along the constrained dynamics trajectory

  • Any qutrit state diagonal in the computational basis

This covers the states encountered in the inaccessible game analysis.

Computational Benefit

With block structure, eigenvalues come from a quadratic (not cubic) formula:

\[\lambda_{1,2} = \frac{(a+c) \pm \sqrt{(a-c)^2 + 4b^2}}{2}, \quad \lambda_3 = d\]

This makes symbolic differentiation ~100× faster than the general case.

References

  • Byrd & Khaneja, Phys. Rev. A 68 (2003)

  • Kimura, Phys. Lett. A 314 (2003)

  • Gamel, Phys. Rev. A 93, 062320 (2016)

Available Methods

Two approaches are available:

  1. LME Exact (qig.symbolic.lme_exact, recommended for LME dynamics):

    • Uses block decomposition: 9×9 → 3×3 + 2×2 + 1×1×4

    • No Taylor approximation - machine precision (~10⁻¹⁵)

    • Works for all 20 block-preserving generators

  2. General su(9) (qig.symbolic.su9_taylor_approximation):

    • Uses Taylor expansion for exp(K)

    • ~1% error at order 2, ~0.0008% at order 6

    • Works for all 80 generators

Usage Example: LME Exact

from qig.symbolic.lme_exact import (
    exact_exp_K_lme,
    exact_rho_lme,
    exact_constraint_lme,
    exact_psi_lme,
    block_preserving_generators,
)
import sympy as sp

# Get available generators
generators, names = block_preserving_generators()
print(f"20 block-preserving generators: {names[:4]}...")

# Create symbolic parameters
a = sp.Symbol('a', real=True)  # local
c = sp.Symbol('c', real=True)  # entangling
theta = {'λ3⊗I': a, 'λ1⊗λ1': c}
theta_list = [a, c]

# EXACT exp(K) - no Taylor approximation!
exp_K = exact_exp_K_lme(theta)

# EXACT constraint C = h₁ + h₂
C = exact_constraint_lme(theta)

# Constraint gradient a = ∇C
a_vec = sp.Matrix([sp.diff(C, t) for t in theta_list])

# Fisher information G = ∇²ψ
psi = exact_psi_lme(theta)
G = sp.Matrix([[sp.diff(sp.diff(psi, ti), tj)
                for tj in theta_list] for ti in theta_list])

# Lagrange multiplier ν = (aᵀGθ)/(aᵀa)
theta_vec = sp.Matrix(theta_list)
nu = (a_vec.T * G * theta_vec)[0,0] / (a_vec.T * a_vec)[0,0]

# Antisymmetric part A = (1/2)[a(∇ν)ᵀ - (∇ν)aᵀ]
grad_nu = sp.Matrix([sp.diff(nu, t) for t in theta_list])
A = (a_vec * grad_nu.T - grad_nu * a_vec.T) / 2

Numeric-Symbolic Bridge

The numeric_lme_blocks_from_theta function bridges the numeric exponential family parameters to the symbolic block structure:

from qig.symbolic.lme_exact import numeric_lme_blocks_from_theta
from qig.exponential_family import QuantumExponentialFamily

# Create numeric exponential family for qutrit pair
qef = QuantumExponentialFamily(n_pairs=1, d=3, pair_basis=True)

# Get Bell state parameters (regularized with log_epsilon for stability)
theta = qef.get_bell_state_parameters(log_epsilon=-20)

# Extract blocks in LME basis
blocks = numeric_lme_blocks_from_theta(theta, qef.operators)

# blocks contains:
# - 'ent_3x3': 3×3 entangled block of K(θ)
# - 'block_2x2': 2×2 block
# - 'diag_1', 'diag_2', 'diag_3', 'diag_4': 1×1 diagonal entries

This enables direct comparison between numeric computations and symbolic expressions for the same state. See the lme_numeric_symbolic_bridge.ipynb notebook for a detailed tutorial.

Precomputed Expressions

For fast evaluation, precomputed symbolic expressions are available in qig/symbolic/precomputed/. These were generated once and saved to Python files:

from qig.symbolic.precomputed.two_param_chain import (
    a, c,  # symbolic parameters
    G, nu, grad_nu,  # intermediate quantities
    M, S, A,  # Jacobian and its decomposition
)

# Evaluate at specific values (e.g., LME scale)
vals = {a: 2.0, c: 2.0}
nu_val = float(nu.subs(vals))
A_num = A.subs(vals)

Caching

Expensive symbolic computations are cached to disk in qig/symbolic/_cache/. The first run may take seconds to minutes; subsequent runs load instantly.

Validation

Run tests with:

pytest tests/test_lme_exact.py -v

Key validations:

  • exp(K): matches scipy.linalg.expm to ~10⁻¹² ✓

  • Constraint gradient a: matches finite difference ✓

  • Fisher info G: matches numerical Hessian ✓

  • ν for local params: equals -1 (structural identity) ✓

  • ν for entangling params: ≠ -1 ✓

  • A for local params: equals 0 ✓

  • A for entangling params: ≠ 0 (Hamiltonian dynamics!) ✓

  • A antisymmetry: A + Aᵀ = 0 ✓

  • S symmetry: S = Sᵀ ✓

  • M decomposition: M = S + A ✓

Current Status

Complete (CIP-0007):

  • EXACT exp(K) for LME dynamics via block decomposition

  • EXACT density matrix ρ and reduced density matrices ρ₁, ρ₂

  • EXACT marginal entropies h₁, h₂

  • EXACT constraint gradient a = ∇(h₁ + h₂)

  • EXACT Fisher information G = ∇²ψ

  • EXACT Lagrange multiplier ν = (aᵀGθ)/(aᵀa)

  • EXACT gradient ∇ν

  • EXACT antisymmetric part A = (1/2)[a(∇ν)ᵀ - (∇ν)aᵀ]

  • EXACT constraint Hessian ∇²C

  • EXACT (∇G)[θ] term

  • EXACT full Jacobian M = -G - (∇G)[θ] + ν∇²C + a(∇ν)ᵀ

  • EXACT symmetric part S = (M + Mᵀ)/2

  • 20 block-preserving generators identified

  • Precomputed expressions for 2-parameter subset

Key results:

  • Local parameters only: ν = -1, ∇ν = 0, A = 0 (structural identity holds)

  • With entangling parameters: ν ≠ -1, ∇ν ≠ 0, A ≠ 0 (Hamiltonian dynamics!)

  • Results verified at LME scale (||θ|| ~ 3): A ≠ 0 for mixed local+entangling params

Planned:

  • Qubit (d=2) implementation

  • Extraction of effective Hamiltonian H_eff from A

  • Extraction of diffusion operator D[ρ] from S

See Also