torch_vn_algebra.algebra¶

Source: torch_vn_algebra/algebra.py. Signatures and docstrings as in the source.

Module documentation¶

Type I von Neumann algebra M = ⊕_{c=1}^{C} M_{n_c}(ℂ) with subspace dimensions k_c.

Storage shape: (batch, C, k_max, k_max)
- batch: parallel Monte Carlo samples
- C: number of factors (channels)
- k_max = max(k_c): maximum subspace dimension across channels

Each channel c has an active block of size k_c × k_c in the top-left corner.
The remaining (k_max - k_c) dimensions are padding.

TypeIAlgebra¶

class TypeIAlgebra(n_factors: List[int], k_factors: List[int], hilbert: Optional[HilbertSpace]=None, batch_size: int=1, complex_valued: bool=True, device: Optional[torch.device]=None, precision: str='single', charges: Optional[List]=None)

Type I von Neumann algebra M = ⊕_{c=1}^{C} M_{n_c}(ℂ).

Operator storage shape: (batch, C, k_max, k_max)

TypeIAlgebra.complex_valued¶

@property
def complex_valued(self) -> bool

TypeIAlgebra.precision¶

@property
def precision(self) -> str

TypeIAlgebra.dtype¶

@property
def dtype(self) -> torch.dtype

TypeIAlgebra.like¶

def like(self, complex_valued: Optional[bool]=None, precision: Optional[str]=None, device=None) -> 'TypeIAlgebra'

The same sector structure with another number field / precision / device.

TypeIAlgebra.bytes_per_operator¶

def bytes_per_operator(self, batch_size: int=1) -> int

Memory of one materialised operator of this algebra (padded layout).

TypeIAlgebra.random_isometry¶

def random_isometry(self, n: int, m: int, batch_size: Optional[int]=None) -> torch.Tensor

Haar-random isometry C^m -> C^n (R^m -> R^n for a real algebra), m <= n: the first m columns of a Haar
unitary, drawn directly as the Q factor of an n x m Ginibre matrix with diag(R) > 0.
Returns a tensor of shape (n, m), or (batch_size, n, m) if batch_size is given.

TypeIAlgebra.random_unitary¶

def random_unitary(self, n: int, measure: str='haar', batch_size: Optional[int]=None) -> torch.Tensor

Random n x n unitary (orthogonal if the algebra is real).

measure:
    'haar'    - Haar measure on U(n) (complex) or O(n) (real)
    'haar_su' - Haar measure on SU(n) (complex) or SO(n) (real)
    'coe'     - Circular Orthogonal Ensemble, S = W^T W (complex only)
    'cse'     - Circular Symplectic Ensemble, S = W^R W with
                W^R = J W^T J^T (complex only, n even)
    'diag'    - diagonal matrix of i.i.d. uniform phases (complex only)

Returns a tensor of shape (n, n), or (batch_size, n, n) if batch_size is given.

TypeIAlgebra.Operator¶

class Operator(algebra: 'TypeIAlgebra', generator: Optional[Callable[[], torch.Tensor]]=None, matrix: Optional[torch.Tensor]=None, is_self_adjoint: Optional[bool]=None, is_normal: Optional[bool]=None, is_positive: Optional[bool]=None, is_invertible: Optional[bool]=None, is_projection: Optional[bool]=None, eigenvalues: Optional[List[torch.Tensor]]=None)

Operator in the von Neumann algebra.

Storage: (batch, C, k_max, k_max)

TypeIAlgebra.Operator.matrix¶

@property
def matrix(self) -> torch.Tensor

TypeIAlgebra.Operator.shape¶

@property
def shape(self) -> Tuple[int, ...]

TypeIAlgebra.Operator.is_self_adjoint¶

@property
def is_self_adjoint(self) -> bool

TypeIAlgebra.Operator.is_self_adjoint¶

def is_self_adjoint(self, value: bool)

TypeIAlgebra.Operator.is_normal¶

@property
def is_normal(self) -> bool

TypeIAlgebra.Operator.is_normal¶

def is_normal(self, value: bool)

TypeIAlgebra.Operator.is_positive¶

@property
def is_positive(self) -> bool

TypeIAlgebra.Operator.is_positive¶

def is_positive(self, value: bool)

TypeIAlgebra.Operator.is_invertible¶

@property
def is_invertible(self) -> bool

TypeIAlgebra.Operator.is_invertible¶

def is_invertible(self, value: bool)

TypeIAlgebra.Operator.is_projection¶

@property
def is_projection(self) -> bool

TypeIAlgebra.Operator.is_projection¶

def is_projection(self, value: bool)

TypeIAlgebra.Operator.lambda_max¶

@property
def lambda_max(self) -> torch.Tensor

TypeIAlgebra.Operator.lambda_min¶

@property
def lambda_min(self) -> torch.Tensor

TypeIAlgebra.Operator.trace¶

@property
def trace(self) -> torch.Tensor

TypeIAlgebra.Operator.trace¶

def trace(self, value: torch.Tensor)

TypeIAlgebra.Operator.trace_norm¶

def trace_norm(self) -> torch.Tensor

||A||_1 = sum of the singular values of the active blocks (for self-adjoint A: sum of |eigenvalues|,
computed with the batched Hermitian eigensolver).

TypeIAlgebra.Operator.frobenius_norm¶

def frobenius_norm(self) -> torch.Tensor

(sum_c ||A_c||_F^2)^{1/2}; the padding is zero.

TypeIAlgebra.Operator.operator_norm¶

def operator_norm(self) -> torch.Tensor

TypeIAlgebra.Operator.inverse¶

def inverse(self, tol: Optional[float]=None) -> 'TypeIAlgebra.Operator'

(Pseudo-)inverse by SVD of every block; singular values below tol are treated as zero (default:
relative tolerance k_c * eps * largest singular value of the block, as numpy.linalg.pinv per block).

TypeIAlgebra.Operator.inv¶

@property
def inv(self) -> 'TypeIAlgebra.Operator'

TypeIAlgebra.Operator.abs¶

def abs(self) -> 'TypeIAlgebra.Operator'

|A| = (A^* A)^{1/2}: for self-adjoint A by the functional calculus, otherwise from the singular
value decomposition (on GPUs for blocks > 32 via the Hermitian eigensolver of [[0, A], [A^*, 0]]);
gradients finite at degenerate singular values.

TypeIAlgebra.Operator.polar¶

def polar(self) -> Tuple['TypeIAlgebra.Operator', 'TypeIAlgebra.Operator']

Polar decomposition A = W |A|, block by block: W unitary for invertible A (a partial isometry on the
kernel of a singular A), |A| = (A^* A)^{1/2}. |A| is differentiable; W carries no gradient.

TypeIAlgebra.Operator.sqrt¶

def sqrt(self) -> 'TypeIAlgebra.Operator'

A^{1/2} of a positive operator (spectral calculus on the active blocks).

TypeIAlgebra.Operator.trace_a_log_a¶

def trace_a_log_a(self, tol: float=1e-12) -> torch.Tensor

TypeIAlgebra.Operator.entropy¶

def entropy(self) -> torch.Tensor

TypeIAlgebra.Operator.michelson_contrast¶

@property
def michelson_contrast(self) -> torch.Tensor

Michelson contrast Delta(A) = (lambda_max - lambda_min) / (lambda_max + lambda_min)
of a positive operator, with the extreme eigenvalues taken over the whole
algebra (all channels). Returns a tensor of shape (batch,); Delta = 0 for A = 0.

TypeIAlgebra.Operator.adjoint¶

def adjoint(self) -> 'TypeIAlgebra.Operator'

TypeIAlgebra.Operator.cast¶

def cast(self, algebra: 'TypeIAlgebra') -> 'TypeIAlgebra.Operator'

This operator as an element of `algebra` (same sectors, other field / precision / device),
e.g. alg.like(precision='double'). Casting a complex operator to a real algebra drops the
imaginary part and raises if it is not negligible.

TypeIAlgebra.Operator.eigh¶

def eigh(self) -> List[Tuple[torch.Tensor, torch.Tensor]]

Spectral decomposition of a self-adjoint operator, one entry per channel:
(eigenvalues (batch, k_c), eigenvectors (batch, k_c, k_c)), A_c = V diag(w) V*.
Channels with k_c = 0 give empty tensors.

TypeIAlgebra.Operator.eigenvalues¶

def eigenvalues(self) -> List[torch.Tensor]

Eigenvalues of a self-adjoint operator, one real tensor (batch, k_c) per channel, ascending.
Cheaper than eigh() (no eigenvectors are formed or stored); differentiable.
Channels with k_c = 0 give empty tensors.

TypeIAlgebra.Operator.cache_spectrum¶

def cache_spectrum(self, enable: bool=True) -> 'TypeIAlgebra.Operator'

Keep the eigendecomposition after its first use, so that further functions of this operator
(apply_function, sqrt, log, power, eigenvalues, propagators for many times, ...) cost matrix products
only. Costs one extra copy of the operator. Also switched on for all operators by
`torch_vn_algebra.caching()`. Returns self. The cache is invalidated by in-place changes of the matrix
(PyTorch version counter), but not by changes through `.data`.

TypeIAlgebra.Operator.clear_cache¶

def clear_cache(self)

TypeIAlgebra.Operator.cache_info¶

def cache_info(self) -> dict

TypeIAlgebra.Operator.apply_function¶

def apply_function(self, f: Callable[[torch.Tensor], torch.Tensor], per_block: bool=False) -> 'TypeIAlgebra.Operator'

f(A) by the spectral theorem, f(A)_c = V_c f(w_c) V_c*, for self-adjoint A (f on the real spectrum)
and for normal A such as unitaries (f on the complex spectrum; f must accept complex tensors).
`f` must act elementwise: it is called once per block size with the eigenvalues of all n sectors of
that size, a tensor (batch, n * k_c), and returns a tensor of the same shape (per-sample parameters
of shape (batch, 1) broadcast; the operator must then have the full batch, e.g. (A * t) instead of a
batch-1 A). With per_block=True f instead receives one block per row, (batch * n, k_c), and may
couple the eigenvalues of a block (e.g. thresholds relative to the block); per-sample parameters do
not broadcast then. A complex-valued result requires a complex algebra. Only the Hermitian
eigensolver is used (for normal A, on H + alpha K with A = H + iK), one batched call per block size.
Gradients with respect to A use the Daleckii-Krein formula (divided differences of f): finite also
at degenerate spectra; for normal non-self-adjoint A they assume f holomorphic near the spectrum.
Gradients with respect to parameters of f (e.g. t in exp(-itw)) are propagated as usual. For
non-normal A use holomorphic_function.

TypeIAlgebra.Operator.holomorphic_function¶

def holomorphic_function(self, f: Callable[[torch.Tensor], torch.Tensor], center=None, radius=None, nodes: int=32, max_nodes: int=2048, rtol: Optional[float]=None, atol: float=0.0) -> 'TypeIAlgebra.Operator'

f(A) for any operator A and f holomorphic near the spectrum (holomorphic functional calculus),
f(A)_c = (2 pi i)^{-1} oint f(z) (z - A_c)^{-1} dz on a circle, by the trapezoidal rule with batched
LU solves; the number of nodes is doubled until the result changes by less than rtol (default 1e-10 in
double, 1e-5 in single precision) relative to its size; ValueError if the estimate stays above
sqrt(rtol). Default circle: just around the spectrum of every block (functional.default_contour);
f must be holomorphic on the closed disk and accept complex tensors. For functions with singularities
(log, powers, resolvents) pass center and radius: numbers, or callables mapping a stack of blocks
(n, k, k) to tensors (n,). The estimated error is stored as result.contour_error after evaluation.
The error is relative to max(|f(X_c)|, atol): when f(X) may vanish (sin at X = 0, z^3 at a
nilpotent block) pass atol > 0, a scale for f near the spectrum; the absolute accuracy accepted where
|f(X_c)| < atol is rtol * atol. Differentiable through the solves. For self-adjoint or normal A, apply_function is exact and cheaper;
for exp, expm is more accurate when the spectrum is spread.

TypeIAlgebra.Operator.expm¶

def expm(self, scale: Union[float, complex]=1.0) -> 'TypeIAlgebra.Operator'

exp(scale * A) for any operator (torch.linalg.matrix_exp on each active block).

TypeIAlgebra.Operator.log¶

def log(self, eps: float=1e-12) -> 'TypeIAlgebra.Operator'

log(A) of a positive operator; eigenvalues are clamped to eps.

TypeIAlgebra.Operator.power¶

def power(self, p: float) -> 'TypeIAlgebra.Operator'

A^p of a positive operator (negative p: on the support only).

TypeIAlgebra.Operator.Tr_blunt¶

def Tr_blunt(self) -> torch.Tensor

sum_c Tr(A_c).

TypeIAlgebra.Operator.Tr_norm¶

def Tr_norm(self) -> torch.Tensor

sum_c Tr(A_c) / k_c (normalised trace of every factor, summed).

TypeIAlgebra.Operator.tau_vN¶

def tau_vN(self) -> torch.Tensor

(1/C) sum_c Tr(A_c) / k_c over the non-empty sectors: the faithful normal tracial state with equal
sector weights.

TypeIAlgebra.operator_from_eigenvalues¶

def operator_from_eigenvalues(self, eigenvalue_sampler: Callable[[int], torch.Tensor], batch_size: int=1, unitary_measure: str='haar', force_self_adjoint: Optional[bool]=None, force_positive: Optional[bool]=None, force_normal: Optional[bool]=None, force_invertible: Optional[bool]=None, force_projection: Optional[bool]=None) -> Operator

TypeIAlgebra.operator¶

def operator(self, matrix: torch.Tensor, validate: bool=False, **properties) -> Operator

Wrap a tensor of shape (batch, C, k_max, k_max) as an operator of this algebra.
Keyword arguments (is_self_adjoint, is_positive, ...) tag known properties. They are trusted (no
synchronisation, no decomposition) unless validate=True, which checks them once; a wrong flag selects a
wrong algorithm (e.g. a non-Hermitian matrix tagged is_self_adjoint is treated as its Hermitian part).

TypeIAlgebra.from_blocks¶

def from_blocks(self, blocks: List[torch.Tensor]) -> Operator

Operator from a list of C blocks of shape (batch, k_c, k_c) (or (k_c, k_c)).

TypeIAlgebra.random_unitary_operator¶

def random_unitary_operator(self, batch_size: int=1, measure: str='haar') -> Operator

Block-diagonal unitary U = (+)_c U_c in the algebra, with independent
U_c drawn from `measure` (see random_unitary) for every channel and sample.

TypeIAlgebra.identity¶

def identity(self, batch_size: int=1) -> Operator

The unit of the algebra. Its matrix is a read-only expanded view of a cached (C, k_max, k_max) tensor
(no memory per sample); clone it before modifying it in place.

TypeIAlgebra.zero¶

def zero(self, batch_size: int=1) -> Operator

TypeIAlgebra.central¶

def central(self, scalars: List[Union[float, complex, torch.Tensor]], batch_size: int=1) -> Operator

Central element sum_c s_c 1_c: one number per sector, or a tensor (batch,) per sector for one value per
sample. Tensors keep their gradients; complex values need a complex algebra.