Convex hull approximation¶
Convex hull wiki-link is a general math concept: given a set of points (could be infinite, but let's assume finite for simplicity) in Euclidean plane $\{a_1,a_2,\cdots,a_m\}\subset \mathbb{R}^n$, the convex hull is defined as
$$ \mathrm{conv}(\{a\}) = \{ x: x=\sum_i p_ia_i, p_i\geq 0,\sum_ip_i=1 \}. $$
For infinite set of points $A$, the definition should be the intersection of all convex sets containing $A$.
In quantum information field, the density matrix set is the convex hull of all pure states
$$ \{\rho\}=\mathrm{conv}\left( \{x\in\mathbb{C}^d: \lVert x \rVert_2=1\} \right). $$
and the SEP set is the convex hull of all product states
$$ \mathrm{SEP}= \mathrm{conv}\left( \{x\otimes y: x\in\mathbb{C}^{d_A},y\in\mathbb{C}^{d_B}, \lVert x \rVert_2=\lVert y \rVert_2=1\} \right)$$
As said in previous tutorials, the SEP set is difficult to describe and in this tutorial we will approximate it by convex hull of finite number of product states.
references
- "Detecting entanglement by pure bosonic extension" arxiv-link
- Upper bounds for relative entropy of entanglement based on active learning doi-link
- A Separability-Entanglement Classifier via Machine Learning doi-link
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
try:
import numqi
except ImportError:
%pip install numqi
import numqi
Parametrization¶
In the CHA method, we are going to parametrize all SEP states according to the definition.
$$ \rho=\sum_{i\in[N]} p_i |\psi_{A,i}\psi_{B,i}\rangle\langle \psi_{A,i}\psi_{B,i}| $$
where
$$ p_i\geq 0, \sum_{i\in[N]} p_i=1, \langle \psi_{A,i}|\psi_{A,i}\rangle=\langle \psi_{B,i}|\psi_{B,i}\rangle=1. $$
Importantly, (TODO-ref) proves that for $d_A\times d_B$ bipartite system, at most $N=d_A^2d_B^2$. Similar discussion is also available on physics-stackexchange What is the minimum number of separable pure states needed to decompose arbitrary separable states?
To make this parametrization, we need to find function mapping the trainable parameters from Euclidean space to the space of density matrices $\rho$ in SEP. As for the positive coefficients $p_i$, we can use the softplus and the normalization function to generate
$$ \mathrm{SoftPlus}(\theta)=\log(1+e^\theta) $$
$$ p_i=\frac{\mathrm{SoftPlus}(\theta_i)}{\sum_j\mathrm{SoftPlus}(\theta_j)} $$
where $\theta$ is a vector of length $N$. As for the pure states $\psi_{A,i}$ and $\psi_{B,i}$, we can use the following parametrization
$$ |\psi_{A}\rangle=\frac{\theta^{(r)}+\mathrm{i}\theta^{(i)}}{\lVert \theta^{(r)}+\mathrm{i}\theta^{(i)}\rVert_2} $$
where the real part $\theta^{(r)}$ and the imaginary part $\theta^{(i)}$ are all vectors of length $d_A$. Similarly, we can parametrize $|\psi_{B}\rangle$ with $\theta^{(r)}$ and $\theta^{(i)}$ of length $d_B$. Put all these together, we have the function mapping:
$$ \rho=f(\theta)=\sum_{i\in[N]} p_i |\psi_{A,i}\psi_{B,i}\rangle\langle \psi_{A,i}\psi_{B,i}| $$
with the number of trainable parameters $N(2d_A+2d_B+1)$.
Let's experiment with the CHA model built in numqi.entangle submodule.
dimA = 3
dimB = 4
model_cha = numqi.entangle.AutodiffCHAREE((dimA, dimB), num_state=2*dimA*dimB, distance_kind='ree')
print('trainable parameters shape')
for key,value in model_cha.named_parameters():
print(f'{key}:, {value.shape}')
print('parametrized separable density matrix:')
model_cha.set_dm_target(numqi.random.rand_density_matrix(dimA*dimB)) #just for test
loss = model_cha()
rho_cha = model_cha.dm_torch.numpy()
print('eigenvalue:', np.round(np.linalg.eigvalsh(rho_cha), 3))
print('trace:', np.trace(rho_cha).real)
trainable parameters shape manifold.manifold_p.theta:, torch.Size([24]) manifold.manifold_psiA.theta:, torch.Size([24, 6]) manifold.manifold_psiB.theta:, torch.Size([24, 8]) parametrized separable density matrix: eigenvalue: [0.018 0.024 0.033 0.042 0.06 0.066 0.075 0.099 0.112 0.121 0.161 0.19 ] trace: 1.0
Above, we print out the shape of parameters $\theta$ and the parametrized separable density matrix $\rho$. We verify that $\rho$ is indeed a density matrix (all eigenvalues are non-negative and the trace is 1).
Loss function¶
Different tasks will use different loss function. For now, let's focus on find the "closest" separable state $\sigma$ for a given density matrix $\rho$.
distance_kind="gellmann" A naive idea is the Hilbert-Schmidt norm
$$ \mathcal{L}(\sigma;\rho)= \lVert \rho-\sigma \rVert_{\mathrm{HS}}=\sum_{ij}|\rho_{ij}-\sigma_{ij}|^2 $$
$$ \min_{\sigma\in\mathrm{CHA}} \mathcal{L}(\sigma;\rho) $$
where we use $\mathrm{CHA}$ to denote the parametrization above and we will see later that CHA is a good approximation of the $\mathrm{SEP}$ set. Moreover, this "distance" is equivalent to the Euclidean distance in the vectorization of density matrix $\vec{\rho}$ (see doc-link for more details). Apparently, when the loss function is optimized to zero, we have $\rho=\sigma$ and $\sigma$ is separable which means that we find a separable decomposition for $\rho$. However, this loss function is not convex (with respect to the parameters) and there is no guarantee that the global minimum can be find. It should be noted that this optimized value is NOT entanglement measure as it could increase when we apply local operations and classical communication (LOCC) on $\rho$ which violate the requirement of the entanglement measure.
distance_kind="ree" Another commonly used loss function is the quantum relative entropy wiki/quantum-relative-entropy
$$ \mathcal{L}(\sigma;\rho)=\mathrm{Tr}[\rho\log\rho-\rho\log\sigma] $$
$$ \min_{\sigma\in\mathrm{CHA}}\mathcal{L}(\sigma;\rho) $$
$\mathcal{L}$ is nonnegative and it's zero if and only if $\rho=\sigma$. However, quantum relative entropy is not a distance function as it's asymmetrical with respect to $\rho$ and $\sigma$. Most importantly, it's jointly convex with respect to both $\rho$ and $\sigma$, so we will see some semi-definite programming (SDP) to solve some objective function in the later tutorials. However, it's not convex with respect to our parametrization $\theta$ and we have to use non-convex optimization (gradient-based optimization) to solve it. The optimized (minimized) value is called relative entropy of entanglement (REE).
Let's verify some analytical results of REE, e.g. Werner state and Isotropic state (TODO-link).
# (d=3) 11 seconds
# (d=4) 25 seconds
dim = 3
alpha_list = np.linspace(0, 1, 100, endpoint=False) # alpha=1 is unstable for matrix logarithm
dm_target_list = [numqi.state.Werner(dim, x) for x in alpha_list]
ree_analytical = np.array([numqi.state.get_Werner_ree(dim, x) for x in alpha_list])
model = numqi.entangle.AutodiffCHAREE((dim, dim), distance_kind='ree')
ree_cha = []
for dm_target_i in tqdm(dm_target_list):
model.set_dm_target(dm_target_i)
tmp0 = numqi.optimize.minimize(model, theta0='uniform', tol=1e-12, num_repeat=1, print_every_round=0).fun
ree_cha.append(tmp0)
ree_cha = np.array(ree_cha)
fig, ax = plt.subplots()
ax.plot(alpha_list, ree_analytical, "x", label="analytical")
ax.plot(alpha_list, ree_cha, label="CHA")
ax.set_xlim(min(alpha_list), max(alpha_list))
ax.set_ylim(1e-13, 1)
ax.set_yscale('log')
ax.set_xlabel(r"$\beta$")
ax.set_ylabel("REE")
ax.set_title(f'Werner({dim})')
ax.legend()
ax.grid()
fig.tight_layout()
dim = 3
alpha_list = np.linspace(0, 1, 100, endpoint=False) # alpha=1 is unstable for matrix logarithm
dm_target_list = [numqi.state.Isotropic(dim, x) for x in alpha_list]
ree_analytical = np.array([numqi.state.get_Isotropic_ree(dim, x) for x in alpha_list])
model = numqi.entangle.AutodiffCHAREE((dim, dim), distance_kind='ree')
ree_cha = [] #about 11 seconds
for dm_target_i in tqdm(dm_target_list):
model.set_dm_target(dm_target_i)
tmp0 = numqi.optimize.minimize(model, theta0='uniform', tol=1e-12, num_repeat=1, print_every_round=0).fun
ree_cha.append(tmp0)
ree_cha = np.array(ree_cha)
fig, ax = plt.subplots()
ax.plot(alpha_list, ree_analytical, "x", label="analytical")
ax.plot(alpha_list, ree_cha, label="CHA")
ax.set_xlim(min(alpha_list), max(alpha_list))
ax.set_ylim(1e-13, 1)
ax.set_yscale('log')
ax.set_xlabel(r"$\beta$")
ax.set_ylabel("REE")
ax.legend()
ax.grid()
ax.set_title(f"Isotropic({dim})")
fig.tight_layout()
