"""Density-matrix helpers used by entanglement and LBO workflows."""
import logging
import numpy as np
from scipy.linalg import eigh as array_eigh
from scipy.sparse import csr_matrix, isspmatrix
logger = logging.getLogger(__name__)
__all__ = [
"diagonalize_density_matrix",
"get_projector_for_efficient_density_matrix",
"flux_block_projector",
]
[docs]
def diagonalize_density_matrix(rho: np.ndarray | csr_matrix):
"""Diagonalize a dense or sparse Hermitian density matrix."""
if isinstance(rho, np.ndarray):
rho_eigvals, rho_eigvecs = array_eigh(rho)
elif isspmatrix(rho):
rho_eigvals, rho_eigvecs = array_eigh(rho.toarray())
else:
raise TypeError("rho must be a NumPy array or SciPy sparse matrix.")
return rho_eigvals, rho_eigvecs
[docs]
def get_projector_for_efficient_density_matrix(
rho_eigvals, rho_eigvecs, threshold: float
):
"""Build a projector from the dominant eigenvectors of a density matrix."""
sorted_indices = np.argsort(rho_eigvals)[::-1]
rho_eigvals = rho_eigvals[sorted_indices]
rho_eigvecs = rho_eigvecs[:, sorted_indices]
p_columns = np.sum(rho_eigvals > threshold)
while p_columns < 2:
threshold /= 10
p_columns = np.sum(rho_eigvals > threshold)
logger.info(
"SIGNIFICANT EIGENVALUES %s with threshold %s", p_columns, threshold
)
proj = np.zeros((rho_eigvals.shape[0], p_columns), dtype=complex)
for jj in range(p_columns):
proj[:, jj] = rho_eigvecs[:, jj]
return proj
[docs]
def flux_block_projector(rho, flux_sig, truncation, flux_tol=1e-9):
"""Diagonalize a density matrix within symmetry blocks and keep dominant modes.
Used by the plaquette-block LBO: group the RDM basis rows by a conserved
"signature" (for lattice gauge theories, the block's boundary electric
flux -- a superselection label, so the RDM is block-diagonal in it),
diagonalize each block separately, and keep the eigenmodes above a
probability threshold. Because every kept eigenvector lives inside a single
signature block, any operator diagonal in the signature (e.g. the boundary
electric-field generators) stays diagonal after projection -- so link
symmetries survive the truncation.
Parameters
----------
rho : numpy.ndarray or scipy.sparse matrix
``(D, D)`` Hermitian reduced density matrix (densified if sparse).
flux_sig : numpy.ndarray
``(D, L)`` real signatures; rows with identical signatures form one
symmetry block. ``L = 0`` (no signature) puts every row in one block.
truncation : float
Keep eigenvectors with eigenvalue (probability) strictly above this.
flux_tol : float, optional
Quantization tolerance when grouping rows by signature.
Returns
-------
proj : numpy.ndarray
``(D, k)`` complex column projector, columns ordered by decreasing
eigenvalue; each column is supported on a single signature block.
info : dict
Diagnostics: ``kept_dim`` (k), ``subsys_dim`` (D), ``n_flux_blocks``,
``kept_eigvals``, ``discarded_weight`` (the truncation error,
``1 - sum(kept)``), and the off-block residual of ``rho``
(``offblock_residual_abs`` / ``_rel``) measuring how well the
signature is conserved -- it must be ~machine-eps for the block-wise
eigendecomposition (and ``discarded_weight``) to be exact.
"""
# Large blocks may arrive as a sparse RDM (sparse-cached psi matrix); all the
# bookkeeping below is dense (D is at most a few thousand), so densify first.
if isspmatrix(rho):
rho = rho.toarray()
rho = np.asarray(rho)
n_states = rho.shape[0] # D: number of basis rows of the density matrix
# --- 1) bucket rows by signature: rows sharing a signature form one block ---
if flux_sig.shape[1]:
# quantize (divide by flux_tol, round) so ~1e-15 noise cannot split a block
keys = np.round(flux_sig / flux_tol).astype(np.int64)
else:
keys = np.zeros((n_states, 1), dtype=np.int64) # no signature -> one block
groups = {} # signature bytes -> indices of its rows
for row in range(n_states):
groups.setdefault(keys[row].tobytes(), []).append(row)
# --- 2) off-block residual: check rho is block-diagonal in the signature ---
block_mask = np.zeros((n_states, n_states), dtype=bool) # True = within-block
for idx in groups.values():
block_mask[np.ix_(idx, idx)] = True # mark this block's diagonal square
off_block = np.array(rho, copy=True) # copy rho, then...
off_block[block_mask] = 0.0 # ...zero the legitimate within-block entries
rho_scale = max(float(np.max(np.abs(rho))), 1e-300) # scale for rel. residual
residual_abs = float(np.max(np.abs(off_block))) if n_states else 0.0
# --- 3) diagonalize each block; keep modes (weights) above truncation ---
cols, eigvals_kept = [], [] # kept full-length eigenvectors + their weights
best_weight, best_column = -1.0, None # dominant BLOCK-LOCAL mode (fallback)
for idx in groups.values():
idx = np.asarray(idx)
# weights w are probabilities (Schmidt weights of the block)
w, v = diagonalize_density_matrix(rho[np.ix_(idx, idx)])
top = int(np.argmax(w)) # this block's dominant mode...
if w[top] > best_weight: # ...tracked so the fallback stays within one block
best_weight = float(w[top])
best_column = np.zeros(n_states, dtype=complex)
best_column[idx] = v[:, top]
keep = np.where(w > truncation)[0]
for k in keep:
column = np.zeros(n_states, dtype=complex) # full length-D column...
column[idx] = v[:, k] # ...holding this block-local mode at its rows
cols.append(column)
eigvals_kept.append(float(w[k]))
if not cols:
# Every weight fell below `truncation` (very high threshold). Keep the
# single most-probable mode -- taken from its OWN block, so the returned
# column still lives in one signature block (the invariant the downstream
# link generators rely on). Do NOT diagonalize the full rho: its top
# eigenvector could straddle blocks when weights are degenerate.
cols.append(best_column)
eigvals_kept.append(best_weight)
# --- 4) stack kept columns into the (D, k) projector, biggest weight first ---
proj = np.column_stack(cols)
eigvals_kept = np.asarray(eigvals_kept)
order = np.argsort(eigvals_kept)[::-1]
proj = proj[:, order]
eigvals_kept = eigvals_kept[order]
# --- 5) diagnostics ---
info = {
"kept_dim": int(proj.shape[1]), # k: columns kept after truncation
"subsys_dim": int(n_states), # D: full density-matrix size
"n_flux_blocks": len(groups), # number of signature (superselection) blocks
"kept_eigvals": eigvals_kept, # retained probabilities (Schmidt weights)
"discarded_weight": float(1.0 - np.sum(eigvals_kept)), # truncation error
"offblock_residual_abs": residual_abs, # worst cross-block entry (abs)
"offblock_residual_rel": residual_abs / rho_scale, # ...relative to |rho|max
}
logger.debug(
"flux_block_projector: D=%d kept=%d blocks=%d off-block(rel)=%.1e disc=%.1e",
info["subsys_dim"],
info["kept_dim"],
info["n_flux_blocks"],
info["offblock_residual_rel"],
info["discarded_weight"],
)
return proj, info