torch_vn_algebra.decompose¶
Source: torch_vn_algebra/decompose.py. Signatures and docstrings as in the source.
Module documentation¶
From symmetries or generators to the algebra: numerical Wedderburn decomposition.
Given matrices :math:`S=\{S_1,\dots,S_J\}\subset M_n(\mathbb C)` -- the elements of a finite group
representation, the generators of the Lie algebra of a compact group representation, or any set of
matrices -- this module computes the von Neumann algebra they generate,
.. math:: \mathcal A = S'' = \text{the unital } *\text{-algebra generated by } S \cup S^*,
its **Wedderburn (structure) decomposition** (lecture notes, Theorem 2.6)
.. math::
V^*\mathcal A V=\bigoplus_{c=1}^{C} M_{k_c}(\mathbb C)\otimes 1_{m_c},\qquad
V^*\mathcal A' V=\bigoplus_{c=1}^{C} 1_{k_c}\otimes M_{m_c}(\mathbb C),\qquad
\sum_c k_c m_c = n,
with a unitary :math:`V`, and returns it as library objects: ``Decomposition.algebra`` is the
:class:`~torch_vn_algebra.TypeIAlgebra` :math:`\bigoplus_c M_{k_c}`, ``Decomposition.commutant_algebra()``
is :math:`\bigoplus_c M_{m_c}`, and ``to_algebra`` / ``from_algebra`` translate between
:math:`n\times n` matrices and operators of these algebras.
import torch_vn_algebra.decompose as dc
G = dc.permutation_representation(3, 2) # S_3 permuting three qubits: (6, 8, 8)
D = dc.from_group(G)
D.k_factors, D.multiplicities # (1, 2), (4, 2): pi(S_3)'' = M_1 (+) M_2
A = D.to_algebra(G) # Operator of D.algebra, batch 6: the blocks rho_c(g)
D.from_algebra(A) # (6, 8, 8): the matrices again
D.verify() # residuals of the decomposition
Mathematics
-----------
*Isotypic decomposition.* :math:`\mathbb C^n` is an :math:`\mathcal A`-module. It splits into
irreducible submodules (minimal :math:`\mathcal A`-invariant subspaces); two of them are *equivalent*
if some invertible (equivalently: some non-zero, by Schur's lemma) intertwiner maps one onto the
other. With :math:`C` equivalence classes of dimensions :math:`k_c` occurring :math:`m_c` times,
:math:`\mathbb C^n\cong\bigoplus_c\mathbb C^{k_c}\otimes\mathbb C^{m_c}` and :math:`\mathcal A` acts
as :math:`\bigoplus_c\rho_c(a)\otimes 1_{m_c}` once the copies of each class are written in
*aligned* bases (bases in which every copy carries the same matrices :math:`\rho_c(a)`).
*Why a random element of the commutant finds the irreducible pieces* (Maehara and Murota 2011,
who apply the method of Murota, Kanno, Kojima and Kojima 2010 to the commutant). In the canonical form a Hermitian :math:`K\in\mathcal A'` is
:math:`\bigoplus_c 1_{k_c}\otimes\kappa_c` with Hermitian :math:`\kappa_c\in M_{m_c}`. For a generic
(e.g. Gaussian) :math:`K` all eigenvalues of all :math:`\kappa_c` are simple and distinct across
:math:`c`, so the eigenspaces of :math:`K` are exactly the subspaces
:math:`\mathbb C^{k_c}\otimes|e\rangle` (:math:`e` an eigenvector of :math:`\kappa_c`): irreducible
submodules, one per copy, of dimension :math:`k_c`. Their computation needs no knowledge of
:math:`V`: an eigenspace of :math:`K` is :math:`\mathcal A`-invariant because :math:`K` commutes with
:math:`\mathcal A`. If eigenvalues of different copies coincide or are closer than the clustering
tolerance, an eigenspace is a sum of several irreducible pieces. This is detected without knowing
the pieces: the compression :math:`Q^*YQ` of another random :math:`Y\in\mathcal A'` to the eigenspace
is a generic element of the commutant of :math:`\mathcal A` restricted there, which is scalar iff the
restriction is irreducible; if it is not scalar, its eigenspaces split the subspace further.
*Why a second random element aligns the copies* (the alignment step of Murota, Kanno, Kojima and
Kojima 2010, Sec. 4.2, here with elements of the commutant). For two irreducible submodules with orthonormal
bases :math:`Q_i,Q_j` and a generic :math:`Y\in\mathcal A'`, the block :math:`T=Q_i^*YQ_j` is an
intertwiner: :math:`T\rho_j(a)=\rho_i(a)T` with :math:`\rho_i(a)=Q_i^*aQ_i`. By Schur's lemma
:math:`T=0` if the submodules are inequivalent and :math:`T=y\,U` with :math:`U` unitary and
:math:`y\neq0` (generically) if they are equivalent. So the non-zero blocks give the equivalence
classes, and replacing :math:`Q_j` by :math:`Q_jU^*` (:math:`U` = polar factor of :math:`T`)
makes copy :math:`j` carry the same matrices as copy :math:`i`. Copy :math:`j` of class :math:`c`,
basis vector :math:`i`, becomes column :math:`\text{off}_c+i\,m_c+j` of :math:`V`.
*The commutant.* A random element of :math:`\mathcal A'` needs :math:`\mathcal A'` itself:
* for **generators** (``wedderburn(..., group=False)``): :math:`\mathcal A'=\{X: S_jX=XS_j,\ S_j^*X=XS_j^*\}`
is the null space of the linear map :math:`X\mapsto([H_j,X])_j` with the Hermitian parts
:math:`H_j` of the :math:`S_j` (same commutant, same generated :math:`*`-algebra). The dense
:math:`(Jn^2)\times n^2` null-space problem costs :math:`O(Jn^6)`. We first diagonalise a random
Hermitian element :math:`R=U\Lambda U^*` of :math:`\mathcal A` (a random combination of the
generators and of products of two of them): every :math:`X\in\mathcal A'` commutes with
:math:`R`, so :math:`U^*XU` is block diagonal with respect to the eigenvalue clusters of
:math:`R`, which leaves :math:`N=\sum_\lambda d_\lambda^2` unknowns instead of :math:`n^2`
(for a generic :math:`R`, :math:`N=\sum_c k_cm_c^2`). The null space of the reduced
:math:`(Jn^2)\times N` system is computed by an incremental QR factorisation and an SVD of the
:math:`N\times N` triangular factor: :math:`O(n^3+Jn^2N^2)` time, :math:`O(n^2N)` memory.
Clusters are formed with a generous tolerance; merging distinct eigenvalues only enlarges
:math:`N` and is always safe, because all constraints are imposed afterwards. Gaps below
:math:`10\,\mathrm{tol}\,\max|\lambda|` never separate clusters: an :math:`X` coupling such
eigenvalues can commute with the generators up to ``tol``, and only the rank decision below may
exclude it (otherwise e.g. :math:`\mathrm{diag}(1,1+10^{-12})` would lose its commutant :math:`M_2`).
* for a **finite group** (``wedderburn(..., group=True)``, all elements :math:`\pi(g)` given): the
twirl :math:`P(X)=|G|^{-1}\sum_g\pi(g)X\pi(g)^*` is the orthogonal projection of :math:`M_n` onto
:math:`\pi(G)'`, so a random element of :math:`\mathcal A'` costs :math:`O(|G|n^3)` and no basis of
:math:`\mathcal A'` is ever formed. :math:`\dim\pi(G)'=|G|^{-1}\sum_g|\chi(g)|^2=\sum_cm_c^2` with
the character :math:`\chi(g)=\operatorname{tr}\pi(g)` is used as an exact check.
*Correctness check.* Every decomposition is verified before it is returned: (i) the generators lie
in :math:`\mathcal B=V(\bigoplus_cM_{k_c}\otimes1_{m_c})V^*` (so :math:`\mathcal A\subseteq\mathcal B`);
(ii) fresh random elements of :math:`\mathcal A'` lie in :math:`\mathcal B'`; (iii)
:math:`\sum_cm_c^2=\dim\mathcal A'`. From (i) :math:`\mathcal A'\supseteq\mathcal B'`, and (iii) gives
equality, hence :math:`\mathcal A=\mathcal A''=\mathcal B''=\mathcal B` (double commutant theorem).
A failed check triggers a new random draw (``max_tries``); if all fail, :class:`DecompositionError`
is raised with the residuals.
Tolerances
----------
With :math:`\varepsilon` the machine epsilon of the working precision (``2.2e-16`` double,
``1.2e-7`` single), all generators normalised to unit Frobenius norm and random commutant elements
to unit Frobenius norm:
* ``tol`` (default :math:`\varepsilon^{1/2}`): singular values of the commutator map below
``tol * max(1, sigma_max)`` span the commutant; Hermitian/anti-Hermitian parts of a generator of
relative size below ``tol`` are dropped as rounding noise (consistently: they could only produce
singular values below ``tol``). ``tol`` decides what counts as a symmetry: input known to accuracy
:math:`\delta` needs ``tol`` :math:`>\delta` (e.g. single-precision data decomposed in double).
Singular values within a factor :math:`\varepsilon^{-1/8}` (90 in double) of the threshold
trigger a :class:`DecompositionWarning` (the commutant dimension may depend on ``tol``).
For generators known to accuracy :math:`\eta` (Hermitian parts, unit Frobenius norm), lecture
notes Theorem 4.25: the computed commutant has the exact dimension and is within
:math:`4\eta/\gamma` of the exact one if :math:`2\eta<\theta<\gamma-2\eta`, with :math:`\gamma`
the smallest non-zero singular value of the exact commutator map; the window is empty if
:math:`\eta\ge\gamma/4`, and :math:`\theta<2\eta` can lose the symmetry.
* noise level :math:`\nu` of the commutant samples: :math:`n\varepsilon` for the twirl or a given
basis; for a null-space basis also :math:`\max(\sigma_{\text{null,max}},\varepsilon\max(1,\sigma_{\max}))
/\sigma_{\text{kept,min}}`, an estimate of the distance of the computed commutant from the exact
one (lecture notes, Lemma 4.24(a)). The decision thresholds below use
:math:`\phi=\min(10^3\nu,\nu^{1/2})`, between the noise and the unit scale of the structure
(Theorem 4.28), and :math:`\theta=\max(\varepsilon^{1/2},\phi)`, i.e. :math:`\varepsilon^{1/2}`
for well-conditioned input (``Decomposition.diagnostics['thresholds']``).
* eigenvalue clustering: a gap smaller than :math:`\varepsilon^{1/4}` times the spread of the
spectrum, or than :math:`\phi\,\|K\|`, does not separate clusters. Merging is safe
(merged clusters are split later); a separating gap :math:`\ge\varepsilon^{1/4}` keeps the
eigenspace error below :math:`\varepsilon^{3/4}`.
* irreducibility: a compressed commutant element :math:`Q^*YQ` with
:math:`\|Q^*YQ-\bar y1\|_F\le\theta` is scalar, so the subspace is irreducible; otherwise it is
split by the eigenspaces of :math:`Q^*YQ` (near-degenerate eigenvalues of the first random element
are resolved this way).
* equivalence: an intertwiner block with :math:`\|Q_i^*YQ_j\|_F>\theta` links two irreducible
subspaces; its singular values must agree to :math:`\max(\varepsilon^{1/4},10\theta)` (a multiple of
a unitary, by Schur's lemma).
* acceptance (``verify``): relative residuals :math:`\le\theta`. Achieved residuals are typically
:math:`10^{-15}`--:math:`10^{-13}` in double precision and :math:`10^{-6}` in single precision.
* a-posteriori certificate (``Decomposition.certificate()``, lecture notes Proposition 4.26): the
backward error :math:`r` (distance of the Hermitian generators from the computed algebra
:math:`\mathcal B`), the gap :math:`\gamma_{\mathcal B}` of the commutator map of their projections onto
:math:`\mathcal B` (exactly computable from the blocks), and the window
:math:`(2r,\gamma_{\mathcal B}-2r)` of rank thresholds certified to give :math:`\dim\mathcal A'=\sum_cm_c^2`:
if it is non-empty, the result is the exact decomposition of generators within :math:`r` of the input.
Limitations
-----------
* The decomposition is over :math:`\mathbb C`. For real generators the result is the decomposition of
the complex algebra they generate (the complexification of the real algebra) and :math:`V` is
complex in general. Over :math:`\mathbb R` the simple summands are of real, complex or quaternionic
type (:math:`M_k(\mathbb R)`, :math:`M_k(\mathbb C)`, :math:`M_k(\mathbb H)`); ``real_structure()``
determines the type of every sector and ``real_basis()`` a real orthogonal basis adapted to it (lecture notes,
Section 4.7). Example: the rotation generator ``[[0, -1], [1, 0]]`` generates the field
:math:`\mathbb C` over :math:`\mathbb R`, but :math:`\mathbb C\oplus\mathbb C` over :math:`\mathbb C`
(``k = (1, 1)``, ``m = (1, 1)``): one complex pair.
* The group path needs *all* group elements of a *unitary* representation (the list may contain
each element the same number of times). A compact group is given by its Lie algebra generators
(``wedderburn(generators)``, e.g. :func:`total_spin`); for a connected group they generate the
same algebra as the group (:math:`e^{tX}` is a limit of polynomials in :math:`X` and
:math:`X=\frac{d}{dt}e^{tX}|_0`). For a disconnected compact group add one element of every
component to the generators.
* Dense linear algebra: the generator path is practical up to :math:`n\approx 100`--:math:`300`
(depending on :math:`N`), the group path up to :math:`|G|n^3\approx10^{11}`.
* :math:`V` is unique up to :math:`V\mapsto V\bigoplus_c(u_c\otimes w_c)` with unitaries
:math:`u_c\in U(k_c)`, :math:`w_c\in U(m_c)`, and a permutation of the sectors; the sectors are
returned sorted by :math:`(k_c,m_c)`.
References
----------
* K. Murota, Y. Kanno, M. Kojima, S. Kojima, *A numerical algorithm for block-diagonal decomposition
of matrix \*-algebras with application to semidefinite programming*, Japan J. Indust. Appl. Math.
27 (2010) 125-160 (the alignment of equivalent copies by a second random element, and the reduction by a
random element of the algebra, used here).
* T. Maehara, K. Murota, *A numerical algorithm for block-diagonal decomposition of matrix
\*-algebras with general irreducible components*, Japan J. Indust. Appl. Math. 27 (2010) 263-293.
* T. Maehara, K. Murota, *Algorithm for error-controlled simultaneous block-diagonalization of
matrices*, SIAM J. Matrix Anal. Appl. 32 (2011) 605-620 (the random-commutant-element algorithm used
here).
* W. Fulton, J. Harris, *Representation Theory: A First Course*, Springer GTM 129 (1991): Schur's
lemma, characters (Sec. 1-2), Schur-Weyl duality (Sec. 6.1).
* Lecture notes, man/lecture_notes/02_von_neumann_algebras.md, Theorem 2.6 and Proposition 2.7.
DecompositionError¶
class DecompositionError()
The randomised decomposition failed its verification in every attempt.
DecompositionWarning¶
class DecompositionWarning()
A numerical rank decision was close to its threshold (ill-conditioned input).
twirl¶
def twirl(group: MatrixLike, X: torch.Tensor) -> torch.Tensor
Group average :math:`P(X)=|G|^{-1}\sum_g\pi(g)X\pi(g)^*`.
For a unitary representation given by all its elements, P is the orthogonal (Hilbert-Schmidt)
projection of :math:`M_n` onto the commutant :math:`\pi(G)'` (lecture notes, Section 2.4).
group: (|G|, n, n) or a list of matrices, all elements of the group.
X: (n, n) or (B, n, n). Returns a tensor of the shape of X (dtype promoted with the group's).
Cost O(|G| B n^3); the sum runs over chunks of group elements, memory O(B n^2 * chunk).
commutant_dimension¶
def commutant_dimension(group: MatrixLike) -> int
:math:`\dim\pi(G)'=|G|^{-1}\sum_g|\operatorname{tr}\pi(g)|^2=\sum_cm_c^2` (orthogonality of
characters), for a unitary representation given by all group elements, (|G|, n, n).
Raises ValueError if the value is not an integer (the matrices are not all elements of a unitary
representation).
commutant_basis¶
def commutant_basis(generators: MatrixLike, tol: Optional[float]=None, *, group: bool=False, method: str='auto', seed: Optional[int]=0, precision: Optional[str]=None, return_info: bool=False)
Orthonormal (Hilbert-Schmidt) basis of the commutant
:math:`S'=\{X\in M_n: XS=SX,\ XS^*=S^*X\ \forall S\in\text{generators}\}`.
generators: (J, n, n) tensor or list of (n, n) matrices (real or complex; not necessarily
Hermitian or *-closed: the adjoints are included).
tol: relative singular-value threshold of the null space (default sqrt(eps) of the precision).
group: True if the matrices are *all* elements of a unitary representation of a finite group;
enables the twirl method.
method: 'nullspace' (null space of X -> ([H_j, X])_j after the reduction by a random element of
the algebra; O(n^3 + J n^2 N^2) with N = sum_c k_c m_c^2 for a generic probe, n^6 at worst),
'twirl' (group only: range of the twirl of d + 8 random matrices, d from the characters;
O(d |G| n^3)), or 'auto': the null space (for a group the twirl is never cheaper: by Burnside
d |G| = sum_g |chi(g)|^2 >= chi(e)^2 = n^2); for a group the null space uses three random elements
of span pi(G) = A as generators and its dimension is checked against the characters (falling back
to the twirl).
seed: seed of the random probe (the span does not depend on it); None for fresh entropy.
precision: 'double' or 'single' (default: from the input dtype).
return_info: also return a dict of diagnostics (method, singular values at the threshold).
Returns a tensor (d, n, n) with tr(B_a^* B_b) = delta_ab (complex dtype), and the info dict if asked.
RealStructure¶
class RealStructure(types, partner, J, M, residuals, tol)
Real structure of a decomposition whose algebra is closed under entrywise complex conjugation (for
example generated by real matrices), as returned by :meth:`Decomposition.real_structure`.
Complex conjugation maps sector c onto a sector ``partner[c]``: :math:`\overline{V_c}=V_{d}(J_c\otimes M_c)`
with :math:`d=` ``partner[c]`` and unitary :math:`J_c\in U(k_c)`, :math:`M_c\in U(m_c)` (lecture notes,
Section 4.7). ``types[c]`` is
* ``'R'`` (real type) if d = c and :math:`J_c\bar J_c=1`: the summand is :math:`M_{k_c}(\mathbb R)`;
* ``'C'`` (complex type) if d != c: sectors c and d form one summand :math:`M_{k_c}(\mathbb C)` of the real
algebra (realified: real dimension 2k_c);
* ``'H'`` (quaternionic type) if d = c and :math:`J_c\bar J_c=-1`: the summand is :math:`M_{k_c/2}(\mathbb H)`
(k_c and m_c are even).
Attributes: ``types``, ``partner``, ``J``, ``M`` (lists over sectors), ``residuals`` (dict: 'leak' =
weight of :math:`\overline{V_c}` outside sector d, 'conjugation' = relative distance of the conjugated
generators from the algebra, 'kron' = distance of the conjugation block from
:math:`J\otimes M`, 'sign' = :math:`\max\|J\bar J-\varepsilon\|`), ``tol`` and ``ok``.
Certificate¶
class Certificate(residuals, gap, window, rank, tol_window, links, commutant_dim, certified, trusted, ok)
A-posteriori certificate of a decomposition (lecture notes, Section 4.9), as returned by
:meth:`Decomposition.certificate`. All norms refer to the Hermitian generators :math:`\tilde H_j` the null space is
computed from: the Hermitian and anti-Hermitian parts of the generators, normalised to unit Frobenius norm.
With :math:`\mathcal B=V(\bigoplus_cM_{k_c}\otimes1_{m_c})V^*`, :math:`E_{\mathcal B}` the trace-preserving
conditional expectation and :math:`H_j=E_{\mathcal B}(\tilde H_j)`:
* ``residuals``: dict with 'unitarity' = :math:`\|V^*V-1\|` (operator norm), 'algebra' and 'commutant' (as in
:meth:`Decomposition.verify`: max relative Frobenius distances of the generators from :math:`\mathcal B` and of
commutant samples from :math:`\mathcal B'`), and 'backward' = :math:`r=(\sum_j\|\tilde H_j-H_j\|_F^2)^{1/2}`;
* ``gap``: :math:`\gamma_{\mathcal B}`, the smallest non-zero singular value of the commutator map
:math:`X\mapsto([H_j,X])_j` of the projected generators (``inf`` if that map is zero), computed exactly from the
blocks of the :math:`H_j` (Proposition 4.26(a));
* ``certified``: :math:`\gamma_{\mathcal B}>2r` (and the residuals pass ``verify``). Then the :math:`H_j`, at
distance :math:`r` from the input, generate exactly :math:`\mathcal B` -- the decomposition is exact for them --
and for every threshold :math:`\theta` in ``window`` :math:`=(2r,\gamma_{\mathcal B}-2r)` the commutant of the
input at tolerance :math:`\theta` has dimension :math:`\sum_cm_c^2` (Proposition 4.26, Theorem 4.25);
* ``rank``: for a null-space commutant, dict with 'threshold' (the cut :math:`\theta`), 'sigma_null_max' and
'sigma_kept_min' (the singular values on both sides of the cut), 'gap' = their ratio and 'margin' =
min(threshold / sigma_null_max, sigma_kept_min / threshold) (the singular values are those of the reduced
null-space problem of :func:`commutant_basis`); ``tol_window`` is ``window`` in units of ``tol`` (divided by
max(1, sigma_max)): the tolerances certified by Proposition 4.26(c) for the full commutator map. None otherwise
(twirl or given commutant);
* ``links``: dict with the link 'threshold' of the grouping step, 'tree_min' (the weakest intertwiner norm used
to align copies) and 'cross_max' (the strongest intertwiner norm between pieces put into different sectors),
relative to the unit-norm random element of the commutant;
* ``trusted``: certified, and for a null-space commutant the cut :math:`\theta` lies in ``window``;
* ``ok``: the result of :meth:`Decomposition.verify`.
Decomposition¶
class Decomposition(V: torch.Tensor, k_factors: Sequence[int], multiplicities: Sequence[int], precision: str, diagnostics: Optional[Dict]=None)
Wedderburn decomposition :math:`V^*\mathcal AV=\bigoplus_cM_{k_c}\otimes1_{m_c}` of a von Neumann
algebra :math:`\mathcal A\subseteq M_n(\mathbb C)`, as returned by :func:`wedderburn`.
Attributes
----------
V : (n, n) unitary. Column ``off_c + i * m_c + j`` (``off_c = sum_{d<c} k_d m_d``) is basis vector
``i < k_c`` (the *algebra index*) of copy ``j < m_c`` (the *multiplicity index*) of sector c, i.e.
:math:`V(\bigoplus_c\mathbb C^{k_c}\otimes\mathbb C^{m_c})` in Kronecker order, so that
``V^* a V = blockdiag_c(kron(a_c, eye(m_c)))`` for a in A and
``V^* x V = blockdiag_c(kron(eye(k_c), x_c))`` for x in A'.
k_factors : (k_c), sizes of the simple summands (dimensions of the irreducible representations).
multiplicities : (m_c).
algebra : TypeIAlgebra(k, k): the abstract algebra (+)_c M_{k_c}, on the device and in the precision
of the decomposition.
n, dims = (dim A, dim A') = (sum k_c^2, sum m_c^2), sector_dims = (k_c m_c).
diagnostics : dict (attempts, refinement rounds, residuals, commutant method, ...).
Sectors are sorted by (k_c, m_c). V is unique up to sector permutations and
V -> V (+)_c (u_c (x) w_c) with unitary u_c in U(k_c), w_c in U(m_c).
Decomposition.dim¶
@property
def dim(self) -> int
dim A = sum_c k_c^2.
Decomposition.commutant_dim¶
@property
def commutant_dim(self) -> int
dim A' = sum_c m_c^2.
Decomposition.dims¶
@property
def dims(self) -> Tuple[int, int]
(dim A, dim A') = (sum_c k_c^2, sum_c m_c^2).
Decomposition.sector_dims¶
@property
def sector_dims(self) -> Tuple[int, ...]
Dimensions k_c m_c of the sector subspaces (ranks of the central projections).
Decomposition.offsets¶
@property
def offsets(self) -> Tuple[int, ...]
First column of every sector in V.
Decomposition.commutant_algebra¶
def commutant_algebra(self) -> TypeIAlgebra
TypeIAlgebra(m, m): the abstract commutant (+)_c M_{m_c}.
Decomposition.to_algebra¶
def to_algebra(self, a: MatrixLike) -> Operator
Blocks of matrices of A: a = V (+)_c (a_c (x) 1_{m_c}) V^* -> Operator with blocks a_c.
a: (n, n) or (batch, n, n). Returns an Operator of ``self.algebra`` with batch size ``batch`` (1
for a single matrix), storage (batch, C, k_max, k_max). Vectorised (two n x n products per
matrix and one gather), no host synchronisation. For a matrix outside A the result is the block
form of the trace-preserving conditional expectation E_A(a) (the copies are averaged:
a_c = m_c^{-1} tr_{m_c}[(V^* a V)_{cc}]).
Decomposition.from_algebra¶
def from_algebra(self, op: Union[Operator, torch.Tensor]) -> torch.Tensor
Inverse of to_algebra: Operator (or tensor (batch, C, k_max, k_max)) -> matrices
V (+)_c (a_c (x) 1_{m_c}) V^*, shape (batch, n, n). Vectorised, no host synchronisation.
Decomposition.to_commutant¶
def to_commutant(self, x: MatrixLike) -> Operator
Blocks of matrices of A': x = V (+)_c (1_{k_c} (x) x_c) V^* -> Operator of
``commutant_algebra()`` with blocks x_c (batch, C, m_max, m_max). Outside A' the result is the
block form of E_{A'}(x) (x_c = k_c^{-1} tr_{k_c}[(V^* x V)_{cc}]).
Decomposition.from_commutant¶
def from_commutant(self, op: Union[Operator, torch.Tensor]) -> torch.Tensor
Inverse of to_commutant: blocks x_c -> V (+)_c (1_{k_c} (x) x_c) V^*, shape (batch, n, n).
Decomposition.conditional_expectation¶
def conditional_expectation(self, a: MatrixLike, onto: str='algebra') -> torch.Tensor
The trace-preserving conditional expectation of M_n onto A (onto='algebra') or A'
(onto='commutant'), i.e. the Hilbert-Schmidt orthogonal projection: (batch, n, n).
For a finite group, E_{A'} is the twirl.
Decomposition.central_projections¶
def central_projections(self) -> torch.Tensor
Minimal central projections P_c = V_c V_c^* of A (and of A'), shape (C, n, n).
Decomposition.basis¶
def basis(self, which: str='algebra') -> torch.Tensor
Hilbert-Schmidt orthonormal basis of A (which='algebra', sum k_c^2 elements
V (E_ii' (x) 1_m) V^* / sqrt(m_c)) or of A' (which='commutant', sum m_c^2 elements
V (1_k (x) E_jj') V^* / sqrt(k_c)); shape (dim, n, n).
Decomposition.real_structure¶
def real_structure(self, tol: Optional[float]=None, seed: Optional[int]=0) -> RealStructure
Types (real, complex, quaternionic) of the simple summands of the *real* algebra, for an algebra
closed under entrywise complex conjugation (generated by real matrices, or by a set closed under
conjugation). See :class:`RealStructure` and lecture notes, Section 4.7.
:math:`\overline{V}=V\,P` with the unitary :math:`P=V^*\overline V`; because :math:`\bar a\in\mathcal A`
for :math:`a\in\mathcal A`, P maps sector c onto one sector d and its block is :math:`J\otimes M` by
Schur's lemma; for d = c, :math:`J\bar J=\varepsilon 1` with :math:`\varepsilon=\pm1`. The type is
decided by the sign of :math:`\operatorname{tr}(J\bar J)/k`, which is :math:`\pm1` up to rounding
(a gap of 2), and partners by the sector carrying the weight of :math:`\overline{V_c}`. Cost:
one n x n product and an SVD of size k_c^2 x m_c^2 per sector.
Raises ValueError if the algebra is not closed under conjugation (the conjugate of a sector is
not, up to ``tol``, inside one sector of the same shape).
Decomposition.real_basis¶
def real_basis(self, tol: Optional[float]=None, seed: Optional[int]=0) -> Tuple[torch.Tensor, List[Dict]]
Real orthogonal matrix O (n x n, real dtype) adapted to the real algebra (lecture notes, Section 4.7):
* real sector c: columns :math:`V_c(W_J\otimes W_M)` with :math:`J=W_JW_J^T`, :math:`M=W_MW_M^T`
(Takagi), so that :math:`O_c^TaO_c=\alpha_c\otimes1_{m_c}` with a *real* :math:`\alpha_c\in M_{k_c}(\mathbb R)`;
* complex pair (c, d): columns :math:`\sqrt2\,\mathrm{Re}\,V_c` then :math:`\sqrt2\,\mathrm{Im}\,V_c`, so that
:math:`O^TaO=\begin{pmatrix}\mathrm{Re}\,x&\mathrm{Im}\,x\\-\mathrm{Im}\,x&\mathrm{Re}\,x\end{pmatrix}` with :math:`x=a_c\otimes1_{m_c}`;
* quaternionic sector c: columns :math:`V_cW` with :math:`J\otimes M=WW^T` (a real basis of the sector).
Returns (O, blocks) with one dict per block: type, sectors, columns (slice). Raises ValueError if
the algebra is not closed under conjugation.
Decomposition.verify¶
def verify(self, tol: Optional[float]=None, generators: Optional[MatrixLike]=None, commutant: Optional[MatrixLike]=None, n_samples: int=2) -> Dict
Residuals of the decomposition (relative Frobenius norms):
* 'unitarity': max |V^* V - 1|;
* 'algebra': max_j ||S_j - E_A(S_j)|| / ||S_j|| over the generators (those given to
:func:`wedderburn`, or ``generators``): zero iff V^* S_j V = (+)_c s_c (x) 1_{m_c};
* 'commutant': the same for elements of A' (``commutant``, or ``n_samples`` fresh random
elements of the commutant used by the decomposition): zero iff V^* X V = (+)_c 1 (x) x_c;
* 'commutant_dim': (dim A' from the commutant basis or the characters, sum_c m_c^2).
Together these certify A = V((+)_c M_{k_c} (x) 1_{m_c})V^* (module docstring). Returns the dict
with 'ok' = all residuals <= tol and the dimensions equal; default tol: the acceptance threshold of the
decomposition, max(sqrt(eps), 1e3 * noise of the commutant samples) (= sqrt(eps) for well-conditioned input).
Decomposition.certificate¶
def certificate(self, generators: Optional[MatrixLike]=None, tol: Optional[float]=None, warn: bool=True) -> Certificate
A-posteriori certificate of the decomposition (lecture notes, Proposition 4.26 and Theorem 4.25): the
residuals, the backward error :math:`r` (distance of the normalised Hermitian generators from the computed
algebra :math:`\mathcal B`), the gap :math:`\gamma_{\mathcal B}` of the commutator map of their projections onto
:math:`\mathcal B`, the window :math:`(2r,\gamma_{\mathcal B}-2r)` of rank thresholds that certify
:math:`\dim\mathcal A'=\sum_cm_c^2`, and the singular values at the cut actually used. See :class:`Certificate`.
generators: the generators (default: those given to :func:`wedderburn`). tol: acceptance tolerance of the
residuals (default: that of :meth:`verify`). warn: issue a :class:`DecompositionWarning` if the result is not
``trusted`` (the gap is not certified, or the cut lies outside the certified window).
Cost: one conditional expectation per generator and one SVD of size :math:`Jk_ck_e\times k_ck_e` per pair of
sectors.
wedderburn¶
def wedderburn(generators: MatrixLike, *, group: bool=False, method: str='random', seed: Optional[int]=0, tol: Optional[float]=None, precision: Optional[str]=None, device=None, commutant: Optional[MatrixLike]=None, max_tries: int=3, max_rounds: int=30) -> Decomposition
Wedderburn decomposition of the von Neumann algebra generated by matrices (module docstring).
generators: (J, n, n) tensor or list of (n, n) matrices. Any matrices: the generated algebra is the
unital *-algebra of S and S^* (S'' by the double commutant theorem), so non-*-closed input is
closed automatically. Lie algebra generators of a connected compact group (Hermitian or
anti-Hermitian) generate the same algebra as the group.
group: True if ``generators`` are *all* elements of a unitary representation of a finite group:
the commutant is then sampled by the twirl (O(|G| n^3) per sample, no commutant basis), and
dim A' is taken from the characters.
method: 'random' (Maehara-Murota with the alignment of Murota-Kanno-Kojima-Kojima; the only method).
seed: seed of the random elements (default 0: reproducible); None for fresh entropy. The sectors and
the subspaces do not depend on it, V only up to the unitary freedom (+)_c u_c (x) w_c.
tol: null-space threshold of the commutant (see commutant_basis); default sqrt(eps).
precision: 'double' (complex128, default) or 'single' (complex64; default for float32/complex64
input; tolerances scale with eps, residuals ~1e-6).
device: device of the computation and of the result (default: that of the input).
commutant: optional precomputed basis (d, n, n) of A' (linearly independent, e.g. from commutant_basis);
d is used as dim A' in the verification.
max_tries: independent random draws before DecompositionError; max_rounds: splitting rounds per draw.
Returns a verified :class:`Decomposition`; ``diagnostics`` records attempts, rounds, residuals and the margins of
the grouping step ('link_threshold', 'link_tree_min', 'link_cross_max'). ``Decomposition.certificate()`` gives an
a-posteriori certificate (lecture notes, Proposition 4.26).
Cost: commutant (see commutant_basis) + O(n^3) per random element + O(J n^3) verification.
from_group¶
def from_group(elements: MatrixLike, **kwargs) -> Decomposition
Decomposition of :math:`\pi(G)''` for a unitary representation of a finite group given by all its
elements (|G|, n, n): ``wedderburn(elements, group=True, **kwargs)``. Then
:math:`V^*\pi(g)V=\bigoplus_c\rho_c(g)\otimes1_{m_c}` with the inequivalent irreducible
representations :math:`\rho_c` (``to_algebra(elements)``) of dimensions k_c and multiplicities m_c
(lecture notes, Proposition 2.7).
from_generators¶
def from_generators(generators: MatrixLike, **kwargs) -> Decomposition
Decomposition of the von Neumann algebra generated by arbitrary matrices (e.g. Lie algebra generators
of a compact group, or generators of a finite group): ``wedderburn(generators, group=False, **kwargs)``.
permutation_representation¶
def permutation_representation(n_sites: int, local_dim: int=2, permutations: Union[str, Sequence[Sequence[int]]]='all', dtype: torch.dtype=torch.complex128, device=None) -> torch.Tensor
The representation of the symmetric group :math:`S_L` permuting the tensor factors of
:math:`(\mathbb C^d)^{\otimes L}`:
:math:`\pi(\sigma)(v_0\otimes\dots\otimes v_{L-1})=v_{\sigma^{-1}(0)}\otimes\dots\otimes v_{\sigma^{-1}(L-1)}`,
i.e. the vector in slot i moves to slot sigma(i); it is a homomorphism,
:math:`\pi(\sigma\tau)=\pi(\sigma)\pi(\tau)` with :math:`(\sigma\tau)(i)=\sigma(\tau(i))`.
n_sites: L; local_dim: d. The basis is the Kronecker one (site 0 most significant, as torch.kron).
permutations: 'all' (the L! elements of S_L, in itertools.permutations order: for from_group),
'adjacent' (the L - 1 transpositions (i, i+1), which generate S_L: for from_generators), or a list
of permutations in one-line notation (p[i] = sigma(i), 0-based).
Returns a real-valued tensor (P, d^L, d^L) of permutation matrices in ``dtype``.
collective_operators¶
def collective_operators(local_ops: MatrixLike, n_sites: int) -> torch.Tensor
Diagonal (collective) action of local operators on :math:`(\mathbb C^d)^{\otimes L}`:
:math:`X\mapsto\sum_{i}1\otimes\dots\otimes X_{(i)}\otimes\dots\otimes1`, the derivative of
:math:`g\mapsto g^{\otimes L}`; for Lie algebra generators X_a of a group acting on C^d these are the
generators of the diagonal action. local_ops: (J, d, d) or list; returns (J, d^L, d^L).
spin_matrices¶
def spin_matrices(spin: float=0.5, dtype: torch.dtype=torch.complex128, device=None) -> torch.Tensor
Spin operators (S_x, S_y, S_z) of spin s in the basis m = s, s-1, ..., -s: (3, 2s+1, 2s+1).
total_spin¶
def total_spin(n_sites: int, spin: float=0.5, dtype: torch.dtype=torch.complex128, device=None) -> torch.Tensor
Total spin :math:`S_a=\sum_i s_a^{(i)}` (a = x, y, z) on :math:`(\mathbb C^{2s+1})^{\otimes L}`:
(3, D, D) with D = (2s+1)^L, the generators of the diagonal SU(2) action. By Schur-Weyl duality
(Fulton-Harris Sec. 6.1) the algebra they generate is the commutant of the permutation algebra:
for s = 1/2, ``from_generators(total_spin(L))`` has k = 2J + 1 (total spin J) and m = dimension of
the S_L irrep of the two-row partition (L/2 + J, L/2 - J).