k-Symmetric density matrix¶
This notebooks demonstrate algorithm in the paper
Tapping into Permutation Symmetry for Improved Detection of k-Symmetric Extensions doi-link
For systems $B_1,B_2,\cdots,B_k$ each of dimension $d$, if the density matirx $\rho\in(\mathcal{H}^d)^{\otimes k}$ is symmetric with respect to the permutation group $S_k$, which means $\rho$ is invariant under any permutation of the subsystems:
$$ P_{B_iB_j} \rho P_{B_iB_j}=\rho $$
above, $P_{B_iB_j}$ is the permutation operator that swaps the $i$-th and $j$-th subsystems.
Then, the density matrix can be written as a direct sum of irreducible representations of the permutation group (TODO-formula):
$$ \rho = \oplus_\lambda \rho_{\lambda} $$
In other words, the density matrix can be block-diagonalized by the permutation group.
TODO
- how to write these basis
- how to find these irrep basis
import time
import numpy as np
import matplotlib.pyplot as plt
try:
import numqi
except ImportError:
%pip install numqi
import numqi
Example: 2 qutrits¶
2 qutrits can be divided into Bosonic subspace (dimension of 6) and Fermionic subspace (dimension of 3). Below, we generate a random density matrix that is symmetric under the permutation group $S_2$.
dim = 3
rho = numqi.random.rand_ABk_density_matrix(dimA=1, dimB=dim, kext=2)
print('density matrix shape:', rho.shape)
density matrix shape: (9, 9)
Let's verify it's symmetric under permutation group $S_2$.
rho_perm = rho.reshape(3, 3, 3, 3).transpose(1,0,3,2).reshape(9,9)
print('error:', np.abs(rho-rho_perm).max())
error: 0.0
Then, let's put the density matrix into the irrep basis.
basis_list = numqi.group.symext.get_sud_symmetric_irrep_basis(dim, kext=2)
basis_list = [y for x in basis_list for y in x] #different young tableaux
for ind0,basis in enumerate(basis_list):
print(f'irrep-basis({ind0}), shape={basis.shape}')
# all irrep-basis make a complete basis (unitary matrix)
unitary = np.concatenate(basis_list, axis=0)
assert np.abs(unitary @ unitary.T.conj() - np.eye(9)).max() < 1e-10
rho_irrep = unitary.conj() @ rho @ unitary.T
rho_irrep[np.abs(rho_irrep)<1e-13] = 0 #remove items with abs value less than 1e-10
fig,ax = plt.subplots()
_ = ax.spy(np.abs(rho_irrep))
irrep-basis(0), shape=(6, 9) irrep-basis(1), shape=(3, 9)
Above, white area means 0 element in the matrix, and nonzero elements are in black. The matrix is block-diagonalized by the permutation group.
Example: 3 qutrits¶
3 qutrits can be divided into Bosonic subspace (dimension of 10), Fermionic basic (dimension of 1), and mixed subspace (dimension of 8) of multiplicity 2.
$$ 27 = 3^3 = 1\times 10 + 2\times 8 + 1\times 1 $$
Let's repeat the same procedure for 3 qutrits.
dim = 3
rho = numqi.random.rand_ABk_density_matrix(dimA=1, dimB=dim, kext=3)
basis_list = numqi.group.symext.get_sud_symmetric_irrep_basis(dim, kext=3)
basis_list = [y for x in basis_list for y in x] #different young tableaux
for ind0,basis in enumerate(basis_list):
print(f'irrep-basis({ind0}), shape={basis.shape}')
unitary = np.concatenate(basis_list, axis=0)
rho_irrep = unitary.conj() @ rho @ unitary.T
rho_irrep[np.abs(rho_irrep)<1e-13] = 0 #remove items with abs value less than 1e-10
fig,ax = plt.subplots()
_ = ax.spy(np.abs(rho_irrep))
irrep-basis(0), shape=(10, 27) irrep-basis(1), shape=(8, 27) irrep-basis(2), shape=(8, 27) irrep-basis(3), shape=(1, 27)
As 2 qutrits example, in the symmetric basis, the density matrix is block-diagonalized by the permutation group $S_3$. What's new is that those two $8\times 8$ blocks are exactly the same, corresoinding to the same Young diagram but different Young tableaux.
block1 = rho_irrep[10:18,10:18]
block2 = rho_irrep[18:26,18:26]
diff = block1 - block2
diff[np.abs(diff)<1e-13] = 0
fig,(ax0,ax1,ax2) = plt.subplots(1,3,figsize=(9,4))
ax0.spy(np.abs(block1))
ax1.spy(np.abs(block2))
ax2.spy(np.abs(diff))
ax0.set_title('1st 8x8 block')
ax1.set_title('2nd 8x8 block')
_ = ax2.set_title('difference')
As shown in the third figure, their difference has only zero elements which means these two $8\times 8$ blocks are exactly the same.
Symmetric extension with irrep¶
When we considering the symmetric extension of the density matrix, we can use the irrep basis to simplify the problem. More importantly, we no need to store the irrep basis explicitly, instead, we only need the following 4-dimensional tensor:
$$ B_{\alpha\beta,rs} = \langle r | \mathrm{Tr}_{B^{k-1}} [|\alpha\rangle\langle \beta|] |s\rangle $$
where $|\alpha\rangle,|\beta\rangle$ are the basis of the irrep, and $|r\rangle,|s\rangle$ are the basis of the subsystem $B$. Below, we show this tensor for 2 qutrits.
coeff_list,multiplicity_list = numqi.group.symext.get_symmetric_extension_irrep_coeff(dim=3, kext=2)
for ind0,coeff in enumerate(coeff_list):
tmp0 = len(np.nonzero(np.abs(coeff)>1e-10)[0])
print(f'irrep({ind0}), shape={coeff.shape}, #nonzero={tmp0}')
numqi.group.symext.print_symmetric_extension_irrep_coeff(coeff)
print()
irrep(0), shape=(6, 6, 3, 3), #nonzero=27 B(3.0,3.0,0.0,0.0)=1.0 B(4.0,4.0,1.0,1.0)=1.0 B(5.0,5.0,2.0,2.0)=1.0 B(0.0,3.0,2.0,0.0)=0.7071067811865475 B(0.0,5.0,0.0,2.0)=0.7071067811865475 B(1.0,3.0,1.0,0.0)=-0.7071067811865475 B(1.0,4.0,0.0,1.0)=0.7071067811865475 B(2.0,4.0,2.0,1.0)=0.7071067811865475 B(2.0,5.0,1.0,2.0)=-0.7071067811865475 B(0.0,0.0,0.0,0.0)=0.4999999999999999 B(0.0,0.0,2.0,2.0)=0.4999999999999999 B(0.0,1.0,2.0,1.0)=-0.4999999999999999 B(0.0,2.0,0.0,1.0)=-0.4999999999999999 B(1.0,1.0,0.0,0.0)=0.4999999999999999 B(1.0,1.0,1.0,1.0)=0.4999999999999999 B(1.0,2.0,0.0,2.0)=0.4999999999999999 B(2.0,2.0,1.0,1.0)=0.4999999999999999 B(2.0,2.0,2.0,2.0)=0.4999999999999999 irrep(1), shape=(3, 3, 3, 3), #nonzero=12 B(0.0,0.0,0.0,0.0)=0.4999999999999999 B(0.0,0.0,2.0,2.0)=0.4999999999999999 B(0.0,1.0,2.0,1.0)=0.4999999999999999 B(0.0,2.0,0.0,1.0)=0.4999999999999999 B(1.0,1.0,0.0,0.0)=0.4999999999999999 B(1.0,1.0,1.0,1.0)=0.4999999999999999 B(1.0,2.0,0.0,2.0)=-0.4999999999999999 B(2.0,2.0,1.0,1.0)=0.4999999999999999 B(2.0,2.0,2.0,2.0)=0.4999999999999999
Time benchmark¶
Making use this tensor, we can reduce the problem to a semidefinite program (SDP) with a small number of variables. As an application, we can solve symmetric extension which is an important criteria to detect entanglement in quantum information theory (see numqi.entangle tutorials for more information). Below, we show the time benchmark for the Werner state.
Werner state boundary $\frac{k+d^2-d}{kd^2+d-1}$, the absolute error is within $10^{-8}$ using MOSEK solver ($10^{-4}$ if SCS solver is used).
# the solver used to build github-pages is SCS, not MOSEK (not free)
case_dict = {
2: [6,8,10,16],
3: [2,3],
}
for dim,kext_list in case_dict.items():
dm_werner = numqi.state.Werner(dim, alpha=1)
dm_norm = numqi.gellmann.dm_to_gellmann_norm(dm_werner)
for kext in kext_list:
_ = numqi.group.symext.get_symmetric_extension_irrep_coeff(dim, kext) #build cache
t0 = time.time()
beta = numqi.entangle.symext.get_ABk_symmetric_extension_boundary(dm_werner, (dim,dim), kext, use_boson=False)
alpha_irrep = (beta/dm_norm)*dim/(beta/dm_norm+dim-1)
alpha_analytical = (kext+dim*dim-dim)/(kext*dim+dim-1)
tmp0 = time.time() - t0
print(f'[d={dim},kext={kext}][{tmp0:.3f}s] alpha={alpha_irrep:.6f}, abs(error)={abs(alpha_analytical-alpha_irrep):.5g}')
[d=2,kext=6][0.021s] alpha=0.615384, abs(error)=2.5785e-07 [d=2,kext=8][0.019s] alpha=0.588235, abs(error)=3.3247e-07 [d=2,kext=10][0.022s] alpha=0.571428, abs(error)=1.0838e-07 [d=2,kext=16][0.038s] alpha=0.545454, abs(error)=4.016e-07 [d=3,kext=2][0.032s] alpha=1.000000, abs(error)=6.8207e-08
[d=3,kext=3][0.087s] alpha=0.818169, abs(error)=1.3297e-05
| $(d,k)$ | QETLAB time (s) | irrep time (s) | $\alpha$ |
|---|---|---|---|
| $(2,6)$ | 0.14 | 0.10 | 0.615 |
| $(2,8)$ | 0.19 | 0.16 | 0.588 |
| $(2,10)$ | 12.60 | 0.16 | 0.571 |
| $(2,16)$ | NA | 0.32 | 0.545 |
| $(2,32)$ | NA | 3.18 | 0.523 |
| $(2,64)$ | NA | 51.96 | 0.512 |
| $(3,3)$ | 0.62 | 0.51 | 0.818 |
| $(3,4)$ | 7.96 | 2.38 | 0.714 |
| $(3,5)$ | NA | 11.56 | 0.647 |
| $(3,6)$ | NA | 55.60 | 0.6 |
- computer spec
- AMD R7-5800H, 16 cpu (hyperthreaded), 16GB memory
- momery is the bottleneck, time may be slightly overestimated
- QETLAB time: kext SDP from QETLAB
- cvx solver for MATLAB
- the matlab script to record time is stored in
project/ws_irrep/demo_qetlab_time.m
- irrep time: SDP with irreducible representation doi-link
- python/cvxpy with MOSEK solver, should be comparable to cvx
Counting parameters¶
dim = 3
kext_list = [3] #[3,4,5,6] #too time-consuming for github-pages
for kext in kext_list:
coeff_list,multiplicity_list = numqi.group.symext.get_symmetric_extension_irrep_coeff(dim, kext) #build cache
dim_list = [x.shape[0] for x in coeff_list]
assert sum(x*y for x,y in zip(dim_list,multiplicity_list))==dim**kext
tmp0 = ' + '.join([f'{x}x{y}' for x,y in zip(dim_list,multiplicity_list)])
num_parameter = sum(x*x for x in dim_list)
print(f'{dim}^{kext} = {tmp0}, #parameter={num_parameter}/{dim**(2*kext)}')
3^3 = 10x1 + 8x2 + 1x1, #parameter=165/729
Many-body systems $(B_1B_2\cdots B_k)$ each of dimension $d$ with Permutation symmetry. $m_i$ is the multiplicity of the irrep $\lambda_i$. $s_i$ is the dimension of the irrep $\lambda_i$. $d^k$ is the dimension of the total Hilbert space. In the following table, the second column is the Young tableau of the irrep, and the third column $\sum_i s_i^2$ counts the total number of parameters in the irrep basis. The last column $d^{2k}$ counts the total number of parameters in the full Hilbert space.
| $(d,k)$ | $$d^k=\sum_i m_i\times s_i$$ | $$\sum_i s_i^2$$ | $$d^{2k}$$ |
|---|---|---|---|
| $(3,3)$ | $$3^3 = 10\times 1 + 8\times 2 + 1\times 1$$ | $165$ | $729$ |
| $(3,4)$ | $$3^4 = 15\times 1 + 15\times 3 + 6\times 2 + 3\times 3$$ | $495$ | $6561$ |
| $(3,5)$ | $$3^5 = 21\times 1 + 24\times 4 + 15\times 5 + 6\times 6 + 3\times 5$$ | $1287$ | $59049$ |
| $(3,6)$ | $$3^6 = 28\times 1 + 35\times 5 + 27\times 9 + 10\times 5 + 10\times 10 + 8\times 16 + 1\times 5$$ | $3003$ | $531441$ |
| $(4,3)$ | $$4^3 = 20\times 1 + 20\times 2 + 4\times 1$$ | $816$ | $4096$ |
| $(4,4)$ | $$4^4 = 35\times 1 + 45\times 3 + 20\times 2 + 15\times 3 + 1\times 1$$ | $3876$ | $65536$ |
| $(4,5)$ | $$4^5 = 56\times 1 + 84\times 4 + 60\times 5 + 36\times 6 + 20\times 5 + 4\times 4$$ | $15504$ | $1048576$ |