English · Español · Français · Deutsch · 中文 · 日本語 · Русский
Tutorial 4. Symmetries and representations¶
This tutorial accompanies Chapter 4 of the lecture notes, which develops the representation
theory behind the algebras of symmetries, turns it into a randomised algorithm that decomposes the von Neumann algebra
generated by any list of matrices (Proposition 4.9), applies it to \(U(1)\) and \(SU(2)\), and sorts the
sectors of an algebra generated by real matrices into real, complex and quaternionic type
(Theorem 4.20), and analyses the decomposition of noisy generators
(Theorem 4.25). We compute characters, run the decomposition on groups and on generators, build spin
chains sector by sector with SpinChain and su2.SU2Chain, compute the real structure of a decomposition, and certify a
decomposition of noisy generators.
4.2 The isotypic decomposition and the commutant¶
Lecture notes: Section 4.2; Theorem 4.3, Theorem 4.4, Corollary 4.5.
Section 4.2 of the lecture notes computes by hand the characters of three irreducible representations of \(S_3\) — the trivial one, the sign and the standard representation on \((1,1,1)^\perp\subseteq\mathbb C^3\) — and of the representation \(\pi\) permuting the three factors of \((\mathbb C^2)^{\otimes3}\), \(\chi_\pi(\sigma)=2^{\#\text{cycles}(\sigma)}\) (table below), and obtains from Corollary 4.5 the multiplicities \(m_{\rm triv}=4\), \(m_{\rm sgn}=0\), \(m_{\rm std}=2\) and \(\dim\pi(S_3)'=20\). The code recomputes the table: it finds \(P_\sigma\) as the restriction of \(\pi(\sigma)\) to the invariant span of \(|100\rangle,|010\rangle,|001\rangle\) (on which \(\pi(\sigma)\) permutes the positions of the single \(1\), Lemma B.1(a)), and \(\#\text{cycles}(\sigma)=\dim\ker(P_\sigma-1)\) (the vectors constant on the cycles).
| identity (\(1\)) | transpositions (\(3\)) | \(3\)-cycles (\(2\)) | |
|---|---|---|---|
| \(\chi_{\rm triv}\) | \(1\) | \(1\) | \(1\) |
| \(\mathrm{sgn}\) | \(1\) | \(-1\) | \(1\) |
| \(\chi_{\rm std}\) | \(2\) | \(0\) | \(-1\) |
| \(\chi_\pi\) | \(8\) | \(4\) | \(2\) |
import math
import torch
import torch_vn_algebra.decompose as dc
torch.manual_seed(0)
dt = torch.complex128
def close(X, Y, tol=1e-10):
return (torch.as_tensor(X) - torch.as_tensor(Y)).abs().max().item() < tol
G = dc.permutation_representation(3, 2) # (6, 8, 8): pi(sigma) permutes three qubits
chi_pi = G.diagonal(dim1=-2, dim2=-1).sum(-1).real # characters Tr pi(sigma)
P3 = G[:, [4, 2, 1]][:, :, [4, 2, 1]].real # on |100>, |010>, |001>: S_3 permuting C^3
cycles = 3 - torch.linalg.matrix_rank(P3 - torch.eye(3, dtype=torch.float64)).double() # dim ker(P - 1)
assert close(chi_pi, 2.0 ** cycles)
irreps = {'trivial': torch.ones(6, dtype=torch.float64), 'sign': (-1.0) ** (3 - cycles),
'standard': P3.diagonal(dim1=-2, dim2=-1).sum(-1) - 1} # fixed points - 1
chis = torch.stack(list(irreps.values()))
assert close(chis @ chis.T / 6, torch.eye(3, dtype=torch.float64)) # Corollary 4.5(a): all irreducible
mult = chis @ chi_pi / 6 # Corollary 4.5(b)
print(dict(zip(irreps, mult.tolist())), (chi_pi ** 2).mean().item())
assert close(mult, torch.tensor([4.0, 0.0, 2.0], dtype=torch.float64))
assert dc.commutant_dimension(G) == 20 # Corollary 4.5(c): M_4 (+) M_2
This is the decomposition of Section 2.4, now obtained from a table of six numbers per representation.
4.3 From generators to the algebra: a randomised algorithm¶
Lecture notes: Section 4.3; Lemma 4.6, Proposition 4.7, Lemma 4.8, Proposition 4.9.
The randomised algorithm of the lecture notes is decompose.wedderburn, with the front ends from_group (all
elements of a finite group; the commutant is sampled with the twirl decompose.twirl) and from_generators (any list
of matrices; the commutant is a null space, decompose.commutant_basis). It splits \(\mathbb C^n\) into the eigenspaces
of a random \(K\in\mathcal A'_h\), groups and aligns them with the blocks of a random \(Y\in\mathcal A'\), and returns a
Decomposition with the unitary V, the \(k_c\) (k_factors) and the \(m_c\) (multiplicities). In floating point it
treats eigenvalues closer than a tolerance as equal and splits the resulting subspaces further by compressing fresh
random elements of the commutant to them (Lemma 4.8(a)).
Decomposition.verify() reports the relative distances of the generators (and of fresh elements of the
commutant) from \(\mathcal B\) (resp. \(\mathcal B'\)), and compares \(\sum_cm_c^2\) with \(\dim\mathcal A'\) (from the
characters for a group, Corollary 4.5(c); from the null space for generators). We run it on \(S_3\) and three
qubits, and then look inside: the eigenspaces of one random \(K\) and the intertwiners of one random \(Y\).
D = dc.from_group(G) # the twirl samples pi(S_3)'
print(D)
assert D.k_factors == (1, 2) and D.multiplicities == (4, 2) # Section 2.4: k = (1, 2), m = (4, 2)
assert D.verify()['ok']
rho = D.to_algebra(G).matrix # (6, 2, 2, 2): the blocks rho_c(sigma), padded
assert close(rho[:, 0, 0, 0], torch.ones(6)) # sector 1: trivial representation
assert close(rho[:, 1].diagonal(dim1=-2, dim2=-1).sum(-1), chis[2]) # sector 2: the standard one
schur = torch.einsum('gij,gkl->ijkl', rho[:, 1], rho[:, 1].conj()) / 6
assert close(schur.reshape(4, 4), torch.eye(4) / 2) # Theorem 4.4 with k = 2
assert close(D.from_algebra(D.to_algebra(G)), G)
Z = torch.randn(8, 8, dtype=dt)
K = dc.twirl(G, Z + Z.mH) # a random element of pi(S_3)'_h
w, Q = torch.linalg.eigh(K)
cuts = [0] + [i + 1 for i in range(7) if w[i + 1] - w[i] > 1e-8] + [8]
spaces = [Q[:, a:b] for a, b in zip(cuts[:-1], cuts[1:])]
assert sorted(W.shape[1] for W in spaces) == [1, 1, 1, 1, 2, 2] # Lemma 4.6: m_c pieces of size k_c
assert all(close(G @ W, W @ (W.mH @ G @ W)) for W in spaces) # invariant subspaces
Y = dc.twirl(G, torch.randn(8, 8, dtype=dt))
Wa, Wb = [W for W in spaces if W.shape[1] == 2]
s = torch.linalg.svdvals(Wa.mH @ Y @ Wb) # Lemma 4.8(b): a multiple of a unitary
assert s[1] > 1e-3 and abs(s[0] - s[1]) < 1e-10
assert all(close(W.mH @ Y @ Wa, 0) for W in spaces if W.shape[1] == 1) # different summands: zero
The commutant dominates the cost: a twirl sample costs \(|G|\) products of \(n\times n\) matrices. For generators
the library first diagonalises a Hermitian element \(R=U\mathrm{diag}(\lambda)U^*\) of \(\mathcal A\) (a random real
combination of the Hermitian parts of the generators plus a multiple of a symmetrised product of two such
combinations). Every \(X\in\mathcal A'\) commutes with \(R\), hence with its spectral projections (Theorem 1.5(d)), so
\(U^*XU\) is block diagonal with respect to the eigenspaces of \(R\), of dimensions \(d_\lambda\); this is imposed in
advance, and the null-space problem has \(N_R=\sum_\lambda d_\lambda^2\) unknowns instead of \(n^2\). Grouping several
eigenvalues together only enlarges the blocks, so the reduced problem has the same solution set \(\mathcal A'\).
In the decomposition, \(V^*RV=\bigoplus_cr_c\otimes1_{m_c}\) with Hermitian \(r_c\in M_{k_c}\); if \(\lambda\) has multiplicity
\(\mu_c(\lambda)\) in \(r_c\), then \(d_\lambda=\sum_c\mu_c(\lambda)m_c\), and since \(\mu^2\ge\mu\) for integers \(\mu\ge0\),
math
N_R=\sum_\lambda\Big(\sum_c\mu_c(\lambda)m_c\Big)^2\ge\sum_{\lambda,c}\mu_c(\lambda)^2m_c^2\ge\sum_{\lambda,c}\mu_c(\lambda)m_c^2=\sum_ck_cm_c^2,
with equality iff the eigenvalues of all the \(r_c\) are simple and pairwise distinct. For the total spin of
\(L=3,4,5\) spins the code below finds equality (numerical observation). Even the minimum is large: for eight spins
\(\frac12\) (\(n=256\)) \(\sum_ck_cm_c^2=4900\), and the reduced system, an \((n^2+N_R)\times N_R\) matrix, has \(3.5\cdot10^8\)
complex entries. For spin chains Section 4.6 therefore uses the module su2, which knows the answer in advance.
for n in (3, 4, 5):
Dn = dc.from_generators(dc.total_spin(n))
_, info = dc.commutant_basis(dc.total_spin(n), return_info=True)
N_min = sum(k * m * m for k, m in zip(Dn.k_factors, Dn.multiplicities))
assert info['unknowns'] >= N_min # N_R >= sum_c k_c m_c^2
print(n, info['unknowns'], N_min) # observed: equality
4.4 \(U(1)\) and charge sectors¶
Lecture notes: Section 4.4; Proposition 4.10.
For \(n\) qubits the charge is the number operator \(\hat N=\sum_jn_j\), \(n_j=|1\rangle\langle1|\) on qubit \(j\), and the operators that conserve it form
\(U(1)'\cong\bigoplus_{w=0}^nM_{\binom nw}(\mathbb C)\) (lecture notes, Section 4.4), which is
chains.SpinChain(n).algebra (bit \(1\) = spin up, \(\hat N=S^z+\frac n2\); in decompose.total_spin spin up is instead
the first basis vector of \(\mathbb C^2\), which relabels the basis but not the sectors). Every XXZ Hamiltonian lies in
\(U(1)'\). The code finds the sectors from \(\hat N\) alone and diagonalises an XXZ ring of four spins sector by sector.
import functools
from torch_vn_algebra import SpinChain
L = 4
Sx, Sy, Sz = dc.total_spin(L) # (16, 16) each, site 0 most significant
N_op = Sz + L / 2 * torch.eye(2 ** L, dtype=dt) # number of up spins
DU = dc.from_generators(N_op) # the algebra generated by U(1)
print(DU)
chain = SpinChain(L, boundary='periodic') # U(1)' = (+)_w M_binom(L, w)
assert set(DU.k_factors) == {1} and sorted(DU.multiplicities) == sorted(chain.algebra.k_factors)
assert sorted(chain.algebra.k_factors) == sorted(math.comb(L, w) for w in range(L + 1))
sx, sy, sz = dc.spin_matrices(0.5)
def two_site(a, b, i, j, n=L): # a on site i, b on site j of n sites
ops = [a if k == i else b if k == j else torch.eye(2, dtype=dt) for k in range(n)]
return functools.reduce(torch.kron, ops)
H = sum(two_site(sx, sx, i, j) + two_site(sy, sy, i, j) + 0.5 * two_site(sz, sz, i, j)
for i, j in chain.bonds) # XXZ ring, Delta = 0.5
assert close(H @ N_op, N_op @ H) # [H, N] = 0 ...
assert close(DU.from_commutant(DU.to_commutant(H)), H) # ... so H is block diagonal (Proposition 4.10)
Hc = chain.xxz(J=1.0, Delta=0.5) # the same H, built sector by sector
E_sec = torch.cat([torch.linalg.eigvalsh(Hc.matrix[0, N, :k, :k]) for N, k in enumerate(chain.algebra.k_factors)])
assert close(torch.sort(E_sec).values, torch.linalg.eigvalsh(H))
4.5 \(SU(2)\): spin¶
Lecture notes: Section 4.5; Theorem 4.11, Proposition 4.13, Theorem 4.14. Spin \(s\) is decompose.spin_matrices(s), and the total spin of \(L\) spins \(\frac12\) is decompose.total_spin(L).
For \(L=4\), lecture notes, Proposition 4.13 gives \(\{S_x,S_y,S_z\}''\cong M_5\oplus M_3\oplus M_1\) with multiplicities \(1,3,2\), of dimension \(35\), and the commutant \(\{S_x,S_y,S_z\}'\cong M_1\oplus M_3\oplus M_2\), of dimension \(14\). The library finds the same from the three generators alone (sectors sorted by \((k_c,m_c)\)), and \(\mathbf S^2=S(S+1)\) on spin \(S\) (Theorem 4.11):
S4 = dc.total_spin(4)
assert close(S4[0] @ S4[1] - S4[1] @ S4[0], 1j * S4[2]) # [S_x, S_y] = i S_z
D4 = dc.from_generators(S4)
print(D4)
assert D4.k_factors == (1, 3, 5) and D4.multiplicities == (2, 3, 1) # M_1, M_3, M_5 with m = 2, 3, 1
assert D4.verify()['ok'] and D4.dims == (35, 14)
m_S = {S: math.comb(4, 2 - S) - (math.comb(4, 1 - S) if S < 2 else 0) for S in (0, 1, 2)}
assert [m_S[(k - 1) // 2] for k in D4.k_factors] == list(D4.multiplicities) # Proposition 4.13
casimir = D4.to_algebra((S4 @ S4).sum(0)) # S^2 = S(S+1) on spin S
for c, k in enumerate(D4.k_factors):
assert close(casimir.matrix[0, c, :k, :k], (k - 1) / 2 * (k + 1) / 2 * torch.eye(k))
By lecture notes, Theorem 4.14 (Schur–Weyl duality) the site permutations span the commutant of the total spin: \(\mathrm{span}\,\pi(S_L)=\{S_x,S_y,S_z\}'\). For \(L\le5\) the code compares the two spaces numerically: by part (a) of the theorem, \(\mathrm{span}\,\pi(S_L)\subseteq\{S_x,S_y,S_z\}'\), it suffices that they have the same dimension.
for n in range(2, 6):
comm = dc.commutant_basis(dc.total_spin(n)).reshape(-1, 4 ** n) # a basis of {S_a}'
perm = dc.permutation_representation(n, 2).reshape(-1, 4 ** n) # all n! permutations
dims = [torch.linalg.matrix_rank(A).item() for A in (comm, perm, torch.cat([comm, perm]))]
print(n, dims)
assert dims[0] == dims[1] == dims[2] # span pi(S_n) = {S_a}'
Dp = dc.from_group(dc.permutation_representation(4, 2)) # pi(S_4)'' = {S_a}'
print(Dp)
assert sorted(zip(Dp.multiplicities, Dp.k_factors)) == sorted(zip(D4.k_factors, D4.multiplicities))
The dimensions printed are \(2,5,14,42\), the Catalan numbers \(\binom{2L}L/(L+1)\); Corollary B.10 proves \(\dim\{S_a\}'=\sum_Sm_S^2=\binom{2L}L/(L+1)\) for every \(L\). The last line is duality in action: \(\pi(S_4)''\) has blocks \(k=(1,2,3)\) with multiplicities \((5,1,3)\) — the pairs \((k_c,m_c)\) of the total spin with the roles exchanged. In general (Corollary B.9) the multiplicity space of spin \(S\) carries an irreducible representation of \(S_L\), of dimension \(m_S\), and these are pairwise inequivalent; here the dimensions are \(1,3,2\) for \(S=2,1,0\).
4.6 Spin systems¶
Lecture notes: Section 4.6; Proposition 4.15, Theorem 4.16, Corollary 4.17.
By lecture notes, Corollary 4.17, orthonormal bases of the highest-weight vectors of \(L\) spins
\(\frac12\) can be built one site at a time with the Clebsch–Gordan coefficients of Theorem 4.16.
Labelling each basis vector by its sequence of intermediate spins \(S_1=\frac12,S_2,\dots,S_L=S\) gives the Yamanouchi
basis that su2.SU2Chain builds (its source uses exactly these coefficients; the code below checks only the
conclusions of Corollary 4.17), inside the \(U(1)\) sector \(\hat N=\frac L2+S\) of SpinChain (whose ordering of the
sites is a relabelling of the tensor factors, which commutes with the total spin by Lemma B.1(b)). The code checks
Theorem 4.16 against dense matrices and against the library's su2.clebsch_gordan, then the conclusions of
Corollary 4.17 for the bases of eight spins. We also check Proposition 4.15(a) on a ring of four spins
with D4, and compare the block sizes for eight spins.
from torch_vn_algebra import su2
def ket(j, m): # |j, m> of Theorem 4.11 (zero if |m| > j)
v = torch.zeros(int(2 * j + 1), dtype=dt)
if abs(m) <= j:
v[int(j - m)] = 1
return v
up, dn = ket(0.5, 0.5), ket(0.5, -0.5)
for j in (0.5, 1.0, 1.5, 2.0): # Theorem 4.16
Sj = dc.spin_matrices(j)
sp_j = [math.sqrt(j * (j + 1) - (j - k) * (j - k + 1)) for k in range(1, int(2 * j) + 1)]
assert close(Sj[0] + 1j * Sj[1], torch.diag(torch.tensor(sp_j, dtype=dt), 1)) # S_+ of Theorem 4.11
assert close(Sj[2], torch.diag(torch.tensor([j - k for k in range(int(2 * j) + 1)], dtype=dt)))
Sc = torch.stack([torch.kron(a, torch.eye(2, dtype=dt)) + torch.kron(torch.eye(len(a), dtype=dt), b)
for a, b in zip(Sj, dc.spin_matrices(0.5))])
Sp, Sm = Sc[0] + 1j * Sc[1], Sc[0] - 1j * Sc[1]
cols = []
for J in (j + 0.5, j - 0.5):
for M in [J - k for k in range(int(2 * J + 1))]:
al, be = math.sqrt((j + M + 0.5) / (2 * j + 1)), math.sqrt((j - M + 0.5) / (2 * j + 1))
p, q = (al, be) if J > j else (-be, al)
x = p * torch.kron(ket(j, M - 0.5), up) + q * torch.kron(ket(j, M + 0.5), dn)
assert close(Sc[2] @ x, M * x)
if M == J:
assert close(Sp @ x, 0) # highest weight
else:
assert close(Sm @ cols[-1], math.sqrt((J + M + 1) * (J - M)) * x) # lowering, (c)
for s2, e in ((0.5, up), (-0.5, dn)): # the library's Racah formula agrees
cg = su2.clebsch_gordan(j, M - s2, 0.5, s2, J, M) if abs(M - s2) <= j else 0.0
assert abs(cg - torch.vdot(torch.kron(ket(j, M - s2), e), x).real.item()) < 1e-12
cols.append(x)
U = torch.stack(cols, 1)
assert U.shape == (2 * len(Sj[0]), 2 * len(Sj[0])) and close(U.mH @ U, torch.eye(U.shape[1])) # (a)
def swap(i, j, n):
p = list(range(n))
p[i], p[j] = j, i
return dc.permutation_representation(n, 2, [p])[0]
ring = [(i, (i + 1) % 4) for i in range(4)]
H4 = sum(swap(i, j, 4) / 2 - torch.eye(16, dtype=dt) / 4 for i, j in ring) # Dirac's identity
assert close(H4, sum(two_site(a, a, i, j) for a in (sx, sy, sz) for i, j in ring))
h = D4.to_commutant(H4) # blocks H_S of sizes m_S = 2, 3, 1
assert close(D4.from_commutant(h), H4) # H4 lies in {S_a}'
levels = torch.cat([torch.linalg.eigvalsh(h.matrix[0, c, :m, :m]).repeat_interleave(k)
for c, (k, m) in enumerate(zip(D4.k_factors, D4.multiplicities))])
assert close(torch.sort(levels).values, torch.linalg.eigvalsh(H4)) # level of H_S repeated 2S+1 times
u1, sym = SpinChain(8, boundary='periodic'), su2.SU2Chain(8, boundary='periodic')
for S, m in zip(sym.spins, sym.multiplicities): # Corollary 4.17: orthonormal highest-weight bases
B, N = sym.hw_basis(S), int(4 + S) # columns in the U(1) sector N = L/2 + S
assert B.shape == (math.comb(8, N), m) and close(B.T @ B, torch.eye(m, dtype=B.dtype))
assert N == 8 or close(sym.raise_(B, N), 0) # S_+ B = 0 (S_z = S holds in sector N)
H_u1, H_su2 = u1.xxz(J=1.0, Delta=1.0), sym.heisenberg(J=1.0)
E_u1 = torch.sort(torch.cat([torch.linalg.eigvalsh(H_u1.matrix[0, N, :k, :k])
for N, k in enumerate(u1.algebra.k_factors)])).values
E_su2, S_lev, _ = sym.full_spectrum(H_su2, return_quantum_numbers=True) # level of block S, 2S+1 times
E_su2, S_lev = E_su2[0], S_lev[0]
assert close(E_u1, E_su2)
assert S_lev[0] == 0 and (S_lev[1:4] == 1).all() and E_su2[4] - E_su2[3] > 1e-3 # singlet, then triplet
assert E_su2[1] - E_su2[0] > 1e-3 and close(E_su2[1:4], E_su2[1].expand(3)) # simple singlet (Theorem 9.9)
assert close(sym.restrict(H_u1).matrix, H_su2.matrix) # Proposition 4.15(b)
flops = lambda ks: sum(k ** 3 for k in ks)
print(u1.algebra.k_factors, sym.multiplicities, flops(u1.algebra.k_factors), flops(sym.multiplicities))
print(f"ground state energy {E_su2[0].item():.6f}, first excited {E_su2[1].item():.6f} (x3)")
For the ring of eight spins the \(U(1)\) blocks have sizes \(\binom8w=1,8,28,56,70,56,28,8,1\) (Section 4.4), the \(SU(2)\) blocks sizes \(m_S=14,28,20,7,1\) for \(S=0,\dots,4\) (Proposition 4.13). The cost of diagonalising, \(\sum k^3\), drops from \(739162\) to \(33040\), by a factor of \(22\); the dimension of the algebra, \(\sum k^2\), from \(\sum_w\binom8w\binom8{8-w}=\binom{16}8=12870\) (choose \(8\) of \(16\) elements, \(w\) of them among the first \(8\)) to \(\sum_Sm_S^2=1430=\binom{16}8/9\) (Corollary B.10). By Theorem 9.9 (the ring has even length and \(J>0\)) the ground state is a single singlet and \(E_0(0)<E_0(1)\). Numerical observation (computed by the code above): the ground energy is \(-3.651093\), and the first excited level is a triplet at \(-3.128419\).
4.7 Real algebras: real, complex and quaternionic type¶
Lecture notes: Section 4.7. For a \(*\)-algebra \(\mathcal A\) closed under complex conjugation, with a decomposition \(V\) whose column blocks \(V_c\) belong to the sectors, Theorem 4.20 gives the permutation \(\sigma\) of the sectors with \(\bar P_c=P_{\sigma(c)}\), unitaries \(J_c\), \(M_c\) with \(\overline{V_c}=V_{\sigma(c)}(J_c\otimes M_c)\) and, for \(\sigma(c)=c\), the sign \(\varepsilon_c\) with \(J_c\bar J_c=\varepsilon_c1\); the types and the real algebra \(\mathcal A_{\mathbb R}\) are in Theorem 4.21, the real Schur lemma is Theorem 4.22, and the real elements of the commutant are described by Proposition 4.23.
Decomposition.real_structure() computes \(\sigma\), \(J_c\), \(M_c\) and the types for a decomposition returned by
wedderburn (from_generators, from_group). It checks that the conjugates of the stored generators lie in
\(V(\bigoplus_cM_{k_c}\otimes1)V^*\) (relative residual 'conjugation') and forms the unitary \(P=V^*\bar V\), whose column block
\(c\) is \(V^*\overline{V_c}\): by Theorem 4.20(a) it is supported on the rows of the sector \(\sigma(c)\), where it equals
\(J_c\otimes M_c\). The library takes for \(\sigma(c)\) the sector carrying the largest part of the block, requires
\((k_{\sigma(c)},m_{\sigma(c)})=(k_c,m_c)\) and bounds the rest ('leak'); it rearranges the block into the
\(k_c^2\times m_c^2\) matrix with entries \(T_{(a,b),(j,l)}=T_{(a,j),(b,l)}\), which for \(T=J\otimes M\) is the rank-one matrix
\(\mathrm{vec}\,J\,(\mathrm{vec}\,M)^T\), reads \(J_c\) and \(M_c\) off its leading singular vectors (with
\(\|J_c\|_F^2=k_c\)), and records the distance of the block from \(J_c\otimes M_c\) ('kron'). If one of these three residuals
exceeds the tolerance, the shapes differ, or the partners do not form an involution, it raises ValueError (the
algebra is not closed under conjugation). For \(\sigma(c)=c\) it sets \(\varepsilon_c\) to the sign of \(\mathrm{Re}\,\mathrm{Tr}(J_c\bar J_c)/k_c\) —
exactly \(\pm1\) in exact arithmetic, so the two candidates are \(2\) apart — and records as
'sign' the maximum over these sectors of \(\|J_c\bar J_c-\varepsilon_c\|_F/\sqrt{k_c}\) and
\(\|M_c\bar M_c-\varepsilon_c\|_F/\sqrt{m_c}\), raised to at least \(1\) if some \(\varepsilon_c=-1\) with odd \(k_c\) or \(m_c\)
(impossible by Theorem 4.20(b)). The result has types, partner (\(\sigma\)),
J, M, residuals and ok, which says that all residuals are at most the tolerance (by default the acceptance
threshold of the decomposition); then \(|\mathrm{Tr}(J_c\bar J_c)/k_c-\varepsilon_c|\le\|J_c\bar J_c-\varepsilon_c\|_F/\sqrt{k_c}\)
(Cauchy–Schwarz) is at most the tolerance too.
Decomposition.real_basis() returns a real orthogonal \(n\times n\) matrix \(O\) and a list that describes its blocks of columns (type, sectors,
column range): for a
sector of real type the columns \(O_c=V_c(W_J\otimes W_M)\) of Theorem 4.21(a), so that \(O_c^TaO_c=\alpha_c(a)\otimes1_{m_c}\)
with real \(\alpha_c(a)\); for a sector of quaternionic type the columns \(V_cW\) with \(J_c\otimes M_c=WW^T\) and
\((J_c\otimes M_c)\bar W=W\) (Lemma 4.19(c) with \(\varepsilon=1\), as \((J_c\otimes M_c)\overline{(J_c\otimes M_c)}=\varepsilon_c^2=1\)),
which are real by the computation of Theorem 4.21(a); and for a pair \(c,\sigma(c)\) (with \(c<\sigma(c)\)) the columns \(\sqrt2\mathrm{Re}\,V_c\) followed by
\(\sqrt2\mathrm{Im}\,V_c\). These are orthonormal: \(V_c^*\overline{V_c}=V_c^*V_{\sigma(c)}(J_c\otimes M_c)=0\) and
\(V_c^T\overline{V_c}=\overline{V_c^*V_c}=1\), with their conjugates, give \((\mathrm{Re}\,V_c)^T\mathrm{Re}\,V_c=(\mathrm{Im}\,V_c)^T\mathrm{Im}\,V_c=\frac12\) and
\((\mathrm{Re}\,V_c)^T\mathrm{Im}\,V_c=0\); columns of different blocks lie in the ranges of the
orthogonal projections \(P_{\mathcal O}=\sum_{e\in\mathcal O}P_e\) of different orbits \(\mathcal O\) of \(\sigma\) (as
\(\overline{V_c}=V_{\sigma(c)}(J_c\otimes M_c)\)), so \(O\) is orthogonal. Taking real and imaginary parts of \(aV_c=V_cx\), \(x=\pi_c(a)\otimes1_{m_c}\),
for real \(a\),
math
O^TaO=\begin{pmatrix}\mathrm{Re}\,x&\mathrm{Im}\,x\\-\mathrm{Im}\,x&\mathrm{Re}\,x\end{pmatrix}\quad\text{on these columns.}
The Takagi factors \(W\) with \(WW^T=J\) of a symmetric unitary \(J=X+iY\) are computed from an orthogonal
diagonalisation of \(X+tY\) for a random \(t\) (\(X\) and \(Y\) are commuting real symmetric matrices, since the imaginary part
of \(J\bar J=JJ^*=1\) is \(YX-XY\)); the library accepts the result on its residuals \(\|WW^T-J\|_F/\|J\|_F\) and
\(\|\mathrm{Im}(V_cW)\|_F/\sqrt{k_cm_c}\) and warns (DecompositionWarning) if one of them exceeds
\(\max(\text{tolerance},10^3\,\epsilon_{\rm mach}\,n)\).
Section 4.7 of the lecture notes works out three algebras generated by real matrices: (i) \(L_a\otimes1_2\) on \(\mathbb R^3\otimes\mathbb R^2\), with the generators \((L_a)_{bc}=-\epsilon_{abc}\) of the rotations of \(\mathbb R^3\): one sector of real type, \(k=3\), \(m=2\); (ii) \(\Omega_2\otimes1_2\) on \(\mathbb R^2\otimes\mathbb R^2\): a pair of complex type, \(k=1\), \(m=2\); (iii) the left multiplications \(1_2\otimes L_q\) by the quaternions on \(\mathbb H^2=\mathbb R^8\): one sector of quaternionic type, \(k=2\), \(m=4\). The code computes the decomposition and the real structure of each, and checks Theorem 4.20(a),(b) and the block form of \(O^TaO\) for the real type (Theorem 4.21(a)) and for the pair.
def qmul(p, q): # H = CD(C), Section E.6: (a, b)(c, d)
(a1, b1), (a2, b2) = p, q
return (a1 * a2 - b2.conjugate() * b1, b2 * a1 + b1 * a2.conjugate())
def q_left(q): # left multiplication by q on H = R^4
cols = [] # coordinates (Re a, Im a, Re b, Im b)
for e in torch.eye(4, dtype=torch.float64).tolist():
u, v = qmul(q, (complex(e[0], e[1]), complex(e[2], e[3])))
cols.append([u.real, u.imag, v.real, v.imag])
return torch.tensor(cols, dtype=torch.float64).T
L_so3 = torch.zeros(3, 3, 3, dtype=torch.float64) # (L_a)_bc = -epsilon_abc
for i, (b1, b2) in enumerate([(1, 2), (2, 0), (0, 1)]):
L_so3[i, b1, b2], L_so3[i, b2, b1] = -1.0, 1.0
R_rot = torch.tensor([[[0.0, -1.0], [1.0, 0.0]]]) # Omega_2
L_H = torch.stack([q_left(u) for u in ((1j, 0), (0, 1), (0, 1j))]) # L_I, L_J, L_K
assert close(L_H[0] @ L_H[1], L_H[2]) and close(q_left((2 - 1j, 3j)).T, q_left((2 + 1j, -3j))) # IJ = K, L_qbar = L_q^T
eye2 = torch.eye(2, dtype=torch.float64)
real_examples = {'R': (torch.kron(L_so3, eye2[None]), (3,), (2,)), # L_a (x) 1_2 on R^3 (x) R^2
'C': (torch.kron(R_rot, eye2[None]), (1, 1), (2, 2)), # Omega_2 (x) 1_2
'H': (torch.kron(eye2[None], L_H), (2,), (4,))} # L_q on H^2 = R^8
dec_real = {}
for t, (gens, ks, ms) in real_examples.items():
D_t = dc.from_generators(gens)
rs_t, (O_t, blocks_t) = D_t.real_structure(), D_t.real_basis()
print(t, D_t, rs_t)
assert D_t.verify()['ok'] and rs_t.ok and (D_t.k_factors, D_t.multiplicities) == (ks, ms)
assert set(rs_t.types) == {t} and close(O_t.T @ O_t, torch.eye(D_t.n)) # real orthogonal basis
for c, d in enumerate(rs_t.partner): # Theorem 4.20
o_c, o_d, w = D_t.offsets[c], D_t.offsets[d], ks[c] * ms[c]
assert close(D_t.V[:, o_c:o_c + w].conj(), D_t.V[:, o_d:o_d + w] @ torch.kron(rs_t.J[c], rs_t.M[c]))
if d == c: # J conj(J) = M conj(M) = +-1
eps = 1.0 if t == 'R' else -1.0
assert close(rs_t.J[c] @ rs_t.J[c].conj(), eps * torch.eye(ks[c]))
assert close(rs_t.M[c] @ rs_t.M[c].conj(), eps * torch.eye(ms[c]))
dec_real[t] = D_t, rs_t, O_t
D_t, rs_t, O_t = dec_real['R']
for g in real_examples['R'][0]: # Theorem 4.21(a): O^T a O = alpha (x) 1_2, alpha real
G_t = O_t.T @ g @ O_t
assert close(G_t, torch.kron(G_t[0::2, 0::2], eye2))
D_t, rs_t, O_t = dec_real['C'] # the pair: [[Re x, Im x], [-Im x, Re x]]
x_t = torch.kron(D_t.to_algebra(real_examples['C'][0]).matrix[0, 0, :1, :1], torch.eye(2, dtype=dt))
assert close(O_t.T @ real_examples['C'][0][0] @ O_t,
torch.cat([torch.cat([x_t.real, x_t.imag], 1), torch.cat([-x_t.imag, x_t.real], 1)]))
The code checks Theorem 4.22 on the irreducible real modules \(\mathbb R^3\), \(\mathbb R^2\), \(\mathbb R^4\) of the three examples (the types and multiplicities found by the library, real with \(m=1\), a complex pair with \(m=1\), quaternionic with \(m=2\), make them irreducible by Theorem 4.22(c)), computing the real commutant as the null space of the real linear equations \(Xs=sX\), \(Xs^T=s^TX\), and then Proposition 4.23 for a random real symmetric element of the commutant in each example.
def real_commutant_dim(gens): # dim_R {X real : X s = s X, X s^T = s^T X}
gens = torch.cat([gens, gens.mT])
eye_n = torch.eye(gens.shape[-1], dtype=torch.float64)
eqs = torch.cat([torch.kron(eye_n, s) - torch.kron(s.T.contiguous(), eye_n) for s in gens])
return gens.shape[-1] ** 2 - torch.linalg.matrix_rank(eqs).item()
for gens, types, ms, dim_D in ((L_so3, ('R',), (1,), 1), (R_rot, ('C', 'C'), (1, 1), 2), (L_H, ('H',), (2,), 4)):
D_irr = dc.from_generators(gens) # Theorem 4.22(c): R^3, R^2, R^4 are irreducible
assert D_irr.real_structure().types == types and D_irr.multiplicities == ms
assert real_commutant_dim(gens) == dim_D # Theorem 4.22(b): D_U = R, C, H
gen_X = torch.Generator().manual_seed(1)
for t, (gens, ks, ms) in real_examples.items(): # Proposition 4.23
D_t, rs_t, O_t = dec_real[t]
B_t = dc.commutant_basis(gens)
X_t = torch.einsum('r,rij->ij', torch.randn(len(B_t), dtype=dt, generator=gen_X), B_t).real
X_t = X_t + X_t.T # real symmetric, in A' (Lemma 4.18(b))
assert close(X_t @ gens, gens @ X_t)
x_t = D_t.to_commutant(X_t).matrix[0]
ev = [torch.linalg.eigvalsh(x_t[c, :m, :m]) for c, m in enumerate(ms)]
if t == 'R': # (a): O^T X O = 1_3 (x) y, y real symmetric
Y_t = O_t.T @ X_t @ O_t
assert close(Y_t, torch.kron(torch.eye(3, dtype=torch.float64), Y_t[:2, :2]))
if t == 'C': # (c): the two sectors have the same spectrum
assert close(ev[0], ev[1])
if t == 'H': # (b): Kramers pairs, each level 2 k_c = 4 times in X
assert close(ev[0][0::2], ev[0][1::2])
assert close(torch.linalg.eigvalsh(X_t), ev[0].repeat_interleave(2))
print('quaternionic block x_c:', [round(v, 4) for v in ev[0].tolist()])
4.8 What to compare with¶
For finite groups, GAP and SageMath compute character tables exactly, and Corollary 4.5 then gives the
multiplicities; RepLAB (MATLAB; D. Rosset, F. Montealegre-Mora, J.-D. Bancal) decomposes representations
of finite and compact groups numerically, also with random elements of the commutant, and
SymbolicWedderburn.jl (Julia) computes Wedderburn decompositions for the symmetry reduction of
semidefinite programs. In physics, QuSpin and other exact-diagonalisation codes block-diagonalise
\(U(1)\) and lattice symmetries, ITensor and TeNPy carry abelian charges in tensor networks, and
TensorKit.jl handles \(SU(2)\) with fusion trees. The library's decompose is a dense method for small
matrices (\(n\) up to a few hundred); it accepts any list of matrices (group elements or generators),
certifies its result (Proposition 4.9), and returns the algebra and its commutant as TypeIAlgebra
objects with to_algebra, from_algebra, to_commutant, from_commutant, so that Chapters 3, 5–8 apply
at once. For spin-\(\frac12\) chains su2.SU2Chain gives the \(SU(2)\) blocks directly (Corollary 4.17); it accepts up to
\(L=18\) spins, a limit set by the memory of the bases.
4.9 Stability of the numerical decomposition¶
Lecture notes: Section 4.9; Lemma 4.24, Theorem 4.25, Proposition 4.26, Proposition 4.27, Theorem 4.28.
The commutant of noisy generators is decided by a tolerance: Theorem 4.25 says
that a threshold \(\theta\) on the singular values of the commutator map \(X\mapsto([H_j,X])_j\) gives the exact commutant
dimension \(d\) if \(2\eta<\theta<\gamma-2\eta\), with \(\eta\) the noise and \(\gamma\) the gap of the exact map, and that the bound
\(2\eta\) on the perturbed zero singular values is attained. The first block computes the singular values of the example of
Theorem 4.25(c), \(\tilde H_1=\mathrm{diag}(1,1,0)+\eta(E_{12}+E_{21})\), and decides the commutant with commutant_basis below,
inside and above the window (the library normalises the generators to unit Frobenius norm, hence the rescaled tol).
import warnings
def commutator_svals(Hs):
"""Singular values, in increasing order, of X -> ([H_j, X])_j with Hilbert-Schmidt norms (row-major vec)."""
n = Hs.shape[-1]
eye = torch.eye(n, dtype=dt)
L = torch.cat([torch.kron(H, eye) - torch.kron(eye, H.T.contiguous()) for H in Hs])
return torch.linalg.svdvals(L).flip(0)
eta = 1e-3
H1 = torch.diag(torch.tensor([1.0, 1.0, 0.0], dtype=dt))
H1t = H1.clone()
H1t[0, 1] = H1t[1, 0] = eta # a symmetry broken at the level eta
s = commutator_svals(H1t[None])
expected = torch.tensor([0, 0, 0, 2 * eta, 2 * eta, 1 - eta, 1 - eta, 1 + eta, 1 + eta], dtype=torch.float64)
assert close(s, expected, 1e-12) # Theorem 4.25(a) and (c): the bound 2 eta is attained
assert close(commutator_svals(H1[None])[4:6], [0, 1]) # exact generator: d = 5, gap 1
scale = torch.linalg.matrix_norm(H1t).item() # theta = tol for unit-norm generators
for theta, dim in ((0.5 * eta, 3), (20 * eta, 5), (2.0, 9)):
with warnings.catch_warnings():
warnings.simplefilter('ignore', dc.DecompositionWarning)
assert dc.commutant_basis(H1t, tol=theta / scale).shape[0] == dim
print('theta < 2 eta: dim 3 (symmetry lost); in the window: 5; above the gap: 9')
Next, the three spin operators of three spins \(\frac12\) (\(k=(2,4)\), \(m=(2,1)\), \(d=5\), gap \(\gamma=1/\sqrt3\) after
normalisation), perturbed by Hermitian noise of increasing size \(\eta\). For each level the code checks Theorem 4.25(a) on
the singular values \(\sigma_d\) and \(\sigma_{d+1}\), puts tol in the middle of the observed window (on a log scale), checks the
dimension and the bound \(\|\tilde\Pi-\Pi\|\le4\eta/\gamma\) of Theorem 4.25(b), and decomposes. With the default
tol \(=\epsilon_{\rm mach}^{1/2}\) the noise is taken seriously as soon as it exceeds the tolerance: the noisy generators generate
\(M_8\), which is the correct answer for them, but not the symmetry one is after. (The DecompositionWarning filtered here
reports that a singular value lies within a factor \(\epsilon_{\rm mach}^{-1/8}\) of the cut.)
S = dc.total_spin(3)
unit = lambda S: S / torch.linalg.matrix_norm(S)[:, None, None] # the Hermitian generators the library uses
sv = commutator_svals(unit(S))
d, gamma = 5, sv[5].item()
assert sv[4] < 1e-12 and abs(gamma - 1 / math.sqrt(3)) < 1e-12
B0 = dc.commutant_basis(S).reshape(d, -1)
Pi0 = B0.T @ B0.conj() # orthogonal projection onto vec(A')
gen = torch.Generator().manual_seed(1)
print(' eta sigma_d sigma_d+1 tol (k, m) default tol')
for eta in (1e-8, 1e-6, 1e-4, 1e-3):
E = torch.randn(S.shape, dtype=dt, generator=gen)
E = E + E.mH
St = S + eta * E / torch.linalg.matrix_norm(E, ord=2)[:, None, None]
eta_n = torch.linalg.matrix_norm(unit(St) - unit(S), ord=2).norm().item() # eta of Theorem 4.25
st = commutator_svals(unit(St))
assert st[d - 1] <= 2 * eta_n and st[d] >= gamma - 2 * eta_n # Theorem 4.25(a)
tol = math.sqrt(st[d - 1] * st[d]) / max(1.0, st[-1].item())
with warnings.catch_warnings():
warnings.simplefilter('ignore', dc.DecompositionWarning)
Bt = dc.commutant_basis(St, tol=tol).reshape(d, -1) # Theorem 4.25(b): dimension d
Dt = dc.from_generators(St, tol=tol)
Dd = dc.from_generators(St)
assert torch.linalg.matrix_norm(Bt.T @ Bt.conj() - Pi0, ord=2) <= 4 * eta_n / gamma
assert (Dt.k_factors, Dt.multiplicities) == ((2, 4), (2, 1))
print(f'{eta_n:.1e} {st[d - 1]:.1e} {st[d]:.3f} {tol:.1e} {Dt.k_factors}, {Dt.multiplicities}'
f' {Dd.k_factors}, {Dd.multiplicities}')
Decomposition.certificate() evaluates Proposition 4.26 on the output: the backward
error \(r\) (the distance of the normalised generators from the computed algebra \(\tilde{\mathcal B}\)), the gap
\(\gamma_{\mathcal B}\) of the commutator map of their projections \(H_j=E(\tilde H_j)\) — computed from the blocks, and checked here
against the dense singular values — and the window \((2r,\gamma_{\mathcal B}-2r)\) of thresholds that certify \(\dim\mathcal A'=\sum_cm_c^2\).
The result is trusted if the threshold actually used lies in that window; otherwise the certificate warns. A tolerance below the
noise gives \(M_8\), certified for these generators, but with a gap of the order of the noise: a model with more symmetry may lie
at distance \(\gamma_{\mathcal B}/2\). A tolerance above the noise but below \(2r\) gives the right structure, which the certificate
cannot vouch for.
c = dc.from_generators(S).certificate()
assert c.certified and c.trusted and c.residuals['backward'] < 1e-12 and abs(c.gap - gamma) < 1e-10
E = torch.randn(S.shape, dtype=dt, generator=gen)
St = S + 1e-4 * (E + E.mH)
with warnings.catch_warnings():
warnings.simplefilter('ignore', dc.DecompositionWarning)
Dt = dc.from_generators(St, tol=1e-2)
D_low = dc.from_generators(St, tol=1e-12) # below the noise
D_mid = dc.from_generators(St, tol=1e-3) # above the noise, below 2 r
ct = Dt.certificate()
print(ct, 'window', [f'{w:.1e}' for w in ct.window], 'threshold', f"{ct.rank['threshold']:.1e}")
assert ct.trusted and 1e-5 < ct.residuals['backward'] < 1e-2 and ct.window[0] < ct.rank['threshold'] < ct.window[1]
sb = commutator_svals(Dt.conditional_expectation(unit(St))) # the projected generators H_j
assert sb[d - 1] < 1e-12 and abs(sb[d] - ct.gap) < 1e-10 # Proposition 4.26(a), (b)
c_low = D_low.certificate()
assert D_low.k_factors == (8,) and c_low.certified and c_low.gap < 1e-3
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter('always')
c_mid = D_mid.certificate()
assert D_mid.k_factors == (2, 4) and not c_mid.trusted and caught
print('tol 1e-12:', D_low.k_factors, f'gap {c_low.gap:.1e}; tol 1e-3:', D_mid.k_factors, 'trusted:', c_mid.trusted)
Finally Proposition 4.27: for a standard Gaussian vector \(G\) of \(\mathcal A'_h\) the blocks are
\(\kappa_c=k_c^{-1/2}H_c\) with independent GUE matrices, and the probability that eigenvalues of different sectors come closer than
\(x\) is at most of order \(x\) (here the bound is \(\sqrt{2/\pi}\,x\,m_1m_2(\frac1{k_1m_1}+\frac1{k_2m_2})^{-1/2}\)), while inside the
\(2\times2\) block the spacing is of order \(x^3\). A Monte Carlo check with \(2\cdot10^4\) samples, drawn in the orthonormal basis
D.basis('commutant') of \(\mathcal A'\):
D = dc.from_generators(S)
coef = torch.randn(20000, d, dtype=dt, generator=gen)
Z = torch.einsum('br,rij->bij', coef, D.basis('commutant')) # standard complex Gaussian vectors of A'
x_c = D.to_commutant((Z + Z.mH) / math.sqrt(2)).matrix # blocks of G: kappa_1 (2 x 2), kappa_2 (1 x 1)
assert abs((x_c[:, 0, :2, :2].abs() ** 2).sum((-2, -1)).mean().item() - 4 / 2) < 0.05 # E Tr kappa_1^2 = m_1^2 / k_1
ev = torch.linalg.eigvalsh(x_c[:, 0, :2, :2])
cross = (ev - x_c[:, 1, :1, 0].real).abs().min(1).values # delta_x
(k1, k2), (m1, m2) = D.k_factors, D.multiplicities
for x in (0.01, 0.03, 0.1):
bound = math.sqrt(2 / math.pi) * x * m1 * m2 / math.sqrt(1 / (k1 * m1) + 1 / (k2 * m2))
p = (cross < x).double().mean().item()
print(f'x = {x}: P(delta_x < x) = {p:.4f} <= {bound:.4f}')
assert p <= bound
p_in = (ev[:, 1] - ev[:, 0] < 0.2).double().mean().item()
assert p_in <= k1 ** 1.5 * 0.2 ** 3 / (6 * math.sqrt(math.pi)) # Proposition 4.27(c), m = 2
Exercises¶
The exercises of the lecture notes, checked numerically (the solutions below are in this order):
1. \(\dim\pi(S_L)'=\binom{d^2+L-1}L\) for \(d=2\) and \(L=3,4\) (the values \(20\) and \(35\));
2. \(1\otimes\frac12\cong\frac12\oplus\frac32\) and \(1\otimes1\cong0\oplus1\oplus2\), with from_generators;
3. the commutant of \(e^{i\theta\hat N}\) on four qubits for \(\theta=1\) (the charges \(0,\dots,4\) give distinct \(e^{iq}\), as
\(0<|q-q'|\le4<2\pi\)) and for \(\theta=\pi\);
4. the Zeeman splitting: the spectrum of \(H+hS_z\) for the Heisenberg ring of eight spins, with su2.SU2Chain and in the
\(U(1)\) sectors;
5. the sectors of the regular representation of \(Q_8\), with real_structure.
Solutions
for n, value in ((3, 20), (4, 35)):
assert dc.commutant_dimension(dc.permutation_representation(n, 2)) == math.comb(n + 3, n) == value
def coupled(a, b): # S^(a) (x) 1 + 1 (x) S^(b)
Sa, Sb = dc.spin_matrices(a), dc.spin_matrices(b)
ea, eb = torch.eye(Sa.shape[-1], dtype=dt), torch.eye(Sb.shape[-1], dtype=dt)
return torch.stack([torch.kron(x, eb) + torch.kron(ea, y) for x, y in zip(Sa, Sb)])
for a, b, ks in ((1.0, 0.5, (2, 4)), (1.0, 1.0, (1, 3, 5))): # 1 (x) 1/2 = 1/2 (+) 3/2, 1 (x) 1 = 0 (+) 1 (+) 2
D_cg = dc.from_generators(coupled(a, b))
assert D_cg.k_factors == ks and set(D_cg.multiplicities) == {1}
for theta, m in ((1.0, (1, 1, 4, 4, 6)), (math.pi, (8, 8))): # theta = 1: e^{iq}, q = 0..4, distinct as 4 < 2 pi
D_theta = dc.from_generators(torch.linalg.matrix_exp(1j * theta * N_op))
assert D_theta.multiplicities == m and set(D_theta.k_factors) == {1}
assert not close(H4 @ Sx - Sx @ H4 + 0.3 * (Sz @ Sx - Sx @ Sz), 0) # [H4 + 0.3 S_z, S_x] != 0
E_field = sym.full_spectrum(sym.heisenberg(J=1.0, h=0.3))[0] # E + h M on every multiplet
H_field = u1.xxz(J=1.0, Delta=1.0, h=0.3)
E_u1_field = torch.cat([torch.linalg.eigvalsh(H_field.matrix[0, N, :k, :k])
for N, k in enumerate(u1.algebra.k_factors)])
assert close(torch.sort(E_u1_field).values, E_field)
Q8 = [(s * a, s * b) for a, b in ((1, 0), (1j, 0), (0, 1), (0, 1j)) for s in (1, -1)]
key = lambda q: tuple(round(v) for z in q for v in (complex(z).real, complex(z).imag))
pos = {key(q): i for i, q in enumerate(Q8)}
reg = torch.zeros(8, 8, 8, dtype=torch.float64) # pi(g) e_h = e_{gh}
for i, g in enumerate(Q8):
for j, h in enumerate(Q8):
reg[i, pos[key(qmul(g, h))], j] = 1.0
D_Q8 = dc.from_group(reg)
rs_Q8 = D_Q8.real_structure()
assert rs_Q8.ok and sorted(zip(rs_Q8.types, D_Q8.k_factors, D_Q8.multiplicities)) == [('H', 2, 2)] + [('R', 1, 1)] * 4
assert D_Q8.dims[0] == 8 # dim_R A_R = dim_C A = |Q_8| (Theorem 4.21)
Next: Tutorial 5: random matrices in an algebra — Haar unitaries, the Weingarten calculus and level statistics in the library.