Geometric measure of subspace entanglement¶
In the previous tutorial (doc-link), we introduce the concept of geometric measure of entanglement $E_G$. In this tutorial, we will introduce how to use numqi to calculate the geometric measure of subspace entanglement. For more details, please check arxiv-link.
import numpy as np
try:
import numqi
except ImportError:
%pip install numqi
import numqi
Recall that the geometric measure of $r$-bounded rank is defined as
$$ E_r(\vert \psi \rangle)=1-\max_{\vert\varphi\rangle \in \sigma_{r-1}}\vert \langle \varphi \vert \psi \rangle\vert^2, $$
where $\sigma_{r-1}$ is the set of all pure states with rank at most $r-1$. Here, rank means the Schmidt rank in the bipartite case, while in the multipartite case, it means the tensor rank. This measure $E_r$ is a generalization of the original one "geometric measure of entanglement" and can be used to quantify different levels of entanglement.
We can go further to define this kind of geometric measure to quantum subspace $\mathcal{S}$. In this case, the geometric measure of subspace entanglement is defined as
$$ E_r(\mathcal{S})=\min_{\vert \psi \rangle \in \mathcal{S}}E_r(\vert \psi \rangle). $$
This measure not only gives a lower bound of the geometric measure for mixed states but also offer a useful approach to investigate the minimal rank of a quantum subspace: if $E_r(\mathcal{S})>0$ and $E_{r+1}(\mathcal{S})=0$, then the minimal rank of $\mathcal{S}$ is $r$. (the definition of minimal rank is given in the previous tutorial doc-link, doc-link).
Method¶
Firstly, we need to modify the definition for further calculation
$$ \begin{aligned} E_r(\mathcal{S}) &= \min _{|\psi\rangle \in \mathcal{S}} E_r(|\psi\rangle) \\ & =\min _{|\psi\rangle \in \mathcal{S}}(1-\max _{\vert \varphi \rangle \in \sigma_{r-1}}|\langle\varphi \vert \psi\rangle|^2) \\ & =1-\max _{\vert \varphi \rangle \in \sigma_{r-1}} \max _{|\psi\rangle \in \mathcal{S}}|\langle\varphi \vert \psi\rangle|^2\\ & =1-\max _{\vert \varphi \rangle \in \sigma_{r-1}}\left\langle\varphi\left|\mathcal{P}_{\mathcal{S}}\right| \varphi\right\rangle \\ & =\min _{\vert \varphi \rangle \in \sigma_{r-1}}\left\langle\varphi\left|\mathcal{P}_{\mathcal{S}}^{\perp}\right| \varphi\right\rangle, \end{aligned} $$ where $\mathcal{P}_{\mathcal{S}}$ projects onto $\mathcal{S}$ and $\mathcal{P}_{\mathcal{S}}^{\perp}$ onto the orthogonal complementary space $\mathcal{S}^{\perp}$. The most crucial step is the transition from the third line to the fourth line follows from the fact that for a given $\vert \varphi \rangle$, the vector from $\mathcal{S}$ maximizing the quantity will be the projection of $\vert \varphi \rangle$ onto $\mathcal{S}$.
Owing to the challenges in characterizing $\sigma_r$ (states with at most tensor rank $r$) and the non-convex nature of optimization, previous work sidestepped a direct computation of $E_r$. Instead, it focused on exploring the lower bound of $E_r$ through mathematical derivation, or by relaxing it in a manner that allows for solution via the semi-definite programming (SDP) method, such as the PPT relaxation:
$$ \begin{aligned} E_G(\mathcal{S})=E_2(\mathcal{S})&=\min _{|\psi_{\text{prod}}\rangle}\left\langle\psi_{\text{prod}}\left|\mathcal{P}_{\mathcal{S}}^{\perp}\right| \psi_{\text{prod}}\right\rangle \\ &=\min _{\left|\psi_{\text {prod}}\right\rangle} \operatorname{tr}\left[\mathcal{P}_{\mathcal{S}}^{\perp}\left|\psi_{\text {prod}}\right\rangle\left\langle\psi_{\text {prod}}\right|\right] \\ & \geqslant \min _{\substack{\rho \geqslant 0 \\ \forall_i \rho^{T_i} \geqslant 0}} \operatorname{tr}\left[\mathcal{P}_{\mathcal{S}}^{\perp} \rho\right], \end{aligned} $$ where $|\psi_{\text{prod}}\rangle$ denotes the \textit{fully} product state in $\sigma_1$ and $\rho^{T_i}$ denotes the partial transpose of $\rho$ with respect to the $i$-th local Hilbert space $\mathcal{H}_i$.
Parametrization and optimization¶
We propose a parametrization strategy to directly characterize the set $\sigma_r$. With the help of parametrization, we can calculate the geometric measure $E_r$ by solving an unconstrained optimization problem. In order to achieve that, we need to find a function $f$ mapping free Euclidean space to the $\sigma_r$.
The states in $\sigma_r$ have the following form:
$$ \vert \psi \rangle = \sum_{i=1}^r \lambda_i \vert \phi_i^{(1)}\rangle \otimes \vert \phi_i^{(2)}\rangle \otimes \cdots \otimes \vert \phi_i^{(n)}\rangle, $$
where the coefficients $\lambda_i$ are positive and $\vert \phi_i^{(k)}\rangle$ is the state in $k$-th local Hilbert space $\mathcal{H}_k$ with dimension $d_k$.
The positive coefficients $\lambda_i$ can be made via a SoftPlus mapping
$$ \lambda=\log\left(1+e^{\theta}\right). $$
As for $\vert \phi_i^{(k)}\rangle$, we can obtain them by the following quotient mapping
$$ |\phi_{i}^{(k)}\rangle=\frac{\alpha_i^{(k)}+\mathrm{i}\beta_{i}^{(k)}}{\left\Vert \alpha_i^{(k)}+\mathrm{i}\beta_{i}^{(k)}\right\Vert _{2}}, $$ where $\alpha_i^{(k)},\beta_i^{(k)}$ are the $k$-th components in $\alpha_i$, $\beta_i$, which are $d_k$-dimensional real vectors.
Example¶
Subspace with analytically known $E_G$¶
The subspace $\mathcal{S}_{2\times d}^{\theta} \subset \mathcal{H}_A \otimes \mathcal{H}_B$ with $d_A=2$ and $d_B=d$, is given by the span of the following vectors
$$ \left|\psi_i\right\rangle_{A B}=a|0\rangle_A\left|i\right\rangle_B+b|1\rangle_A\left|i+1\right\rangle_B, $$ where $i=0,1,\cdots,d-2$, with $a=\cos(\theta/2)$ and $b=\exp(i\xi)\sin(\theta/2)$, $\theta \in (0,\pi)$, $\xi \in [0,2\pi)$.
Clearly, the dimension of $\mathcal{S}_{2\times d}^{\theta}$ is $d-1$, which is also the maximal available dimension of an entangled subspace in this scenario. Furthermore, it has been proved that
$$ E_{2}\left(\mathcal{S}_{2 \times d}^\theta\right)=\frac{1}{2}\left[1-\sqrt{1-\sin ^2 \theta \sin ^2\left(\frac{\pi}{d}\right)}\right]. $$
dimB = 3
theta = np.pi/3
matrix_subspace = numqi.matrix_space.get_GES_Maciej2019(dimB, theta=theta)
# analytical solution
ret_analytical = numqi.matrix_space.get_GM_Maciej2019(dimB, theta)
# PPT aproximation
ret_ppt = numqi.matrix_space.get_generalized_geometric_measure_ppt(matrix_subspace, [2,dimB])
# gradient descent
gd_kwargs = dict(theta0='uniform', num_repeat=3, tol=1e-12, print_every_round=0)
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel((2,dimB), rank=1)
model.set_target(matrix_subspace)
ret_gd = numqi.optimize.minimize(model, **gd_kwargs).fun
print('Analytical solution:', ret_analytical)
print('PPT approximation:', ret_ppt)
print('Gradient descent:', ret_gd)
Analytical solution: 0.16928108611692616 PPT approximation: 0.16928065317845392 Gradient descent: 0.16928108611692594
In the above example, we can see both gradient descent strategy and PPT approximation give the accurate result, consistent with the analytical solution.
Bound entanglement¶
One famous approach to construct bound entangled states is using the unextendible product basis (UPB), as follows
$$ \rho = \frac{1}{d_Ad_B-d_{\mathcal{S}}}\mathcal{P}_{\mathcal{S}}^{\bot}, $$ where $d_A, d_B$ are the local dimensions for bipartite states, $\mathcal{S}$ represents the subspace spanned by the given UPB and $d_{\mathcal{S}}$ is the dimension of that subspace. For example, we can consider the 5-state "Tiles" UPB \cite{divincenzo2003unextendible} (here we omit normalization for brevity),
$$ \begin{aligned} \mathcal{S}_{\text {tiles }}= & \operatorname{span}\{|0\rangle \otimes(|0\rangle-|1\rangle),|2\rangle \otimes(|1\rangle-|2\rangle),\\ &(|0\rangle-|1\rangle) \otimes|2\rangle,(|1\rangle-|2\rangle) \otimes|0\rangle,(|0\rangle+|1\rangle+|2\rangle) \\ & \otimes(|0\rangle+|1\rangle+|2\rangle)\} \subset \mathcal{H}_A \otimes \mathcal{H}_B . \end{aligned} $$ As we have mentioned before, the geometric measure of the support space can give a lower bound for the geometric measure of the state, i.e., $E_r(\rho) \geq E_r(\text{supp}(\rho))$. For the state $\rho_{\text{tiles}}$ constructed by the "Tiles" UPB, we can estimate $E_2(\text{supp}(\rho_{\text{tiles}}))$ by the PPT relaxation and gradient descent and compare the difference between them.
dimA = 3
dimB = 3
rho_bes = numqi.entangle.load_upb('tiles', return_bes=True)[1]
EVL,EVC = np.linalg.eigh(rho_bes)
matrix_subspace = (EVC[:, EVL>1e-4]).reshape(dimA,dimB,-1).transpose(2,0,1)
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel((dimA, dimB), rank=1)
model.set_target(matrix_subspace)
kwargs = dict(theta0='uniform', num_repeat=3, tol=1e-12, print_every_round=0)
theta_optim = numqi.optimize.minimize(model,**kwargs)
print(f'lower bound given by gradient descent = {theta_optim.fun}')
print(f'lower bound given by PPT = {numqi.matrix_space.get_geometric_measure_ppt(matrix_subspace, [dimA, dimB])}')
lower bound given by gradient descent = 0.028416213335729745 lower bound given by PPT = 6.016892294009278e-11
We can see that $E_2(\text{supp}(\rho_{\text{tiles}}))$ calculated by the PPT relaxation is around $10^{-12}$ close to $0$ while the gradient descent gives around $0.0284$. That means PPT cannot detect the entanglement in this case while the gradient descent does.
Higher entangled state¶
Let $d_A=d_B=4$ and consider the following mixed state from
$$ \rho=\frac{1}{3} \sum_{i=1}^3\left|\psi_i\right\rangle\left\langle \psi_i\right| \in \mathcal{H}_A \otimes \mathcal{H}_B, $$ where (we omit normalization for brevity) $$ \begin{aligned} & \left|\psi_1\right\rangle=|0\rangle \otimes|0\rangle+|1\rangle \otimes|1\rangle+|2\rangle \otimes|2\rangle+|3\rangle \otimes|3\rangle, \\ & \left|\psi_2\right\rangle=|0\rangle \otimes|1\rangle+|1\rangle \otimes|2\rangle+|2\rangle \otimes|3\rangle+|3\rangle \otimes|0\rangle, \\ & \left|\psi_3\right\rangle=|0\rangle \otimes|2\rangle+|1\rangle \otimes|3\rangle+|2\rangle \otimes|0\rangle-|3\rangle \otimes|1\rangle . \end{aligned} $$
Similarly, let us calculate the lower bound of $E_2(\rho)$ by PPT relaxation and gradient descent.
dimA = 4
dimB = 4
tmp0 = [
[(0,0,1), (1,1,1), (2,2,1), (3,3,1)],
[(0,1,1), (1,2,1), (2,3,1), (3,0,1)],
[(0,2,1), (1,3,1), (2,0,1), (3,1,-1)],
]
matrix_subspace = np.stack([numqi.matrix_space.build_matrix_with_index_value(dimA, dimB, x) for x in tmp0])
print(f'lower bound given by PPT = {numqi.matrix_space.get_geometric_measure_ppt(matrix_subspace, [dimA, dimB])}')
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel((dimA, dimB), rank=1)
model.set_target(matrix_subspace)
kwargs = dict(theta0='uniform', num_repeat=3, tol=1e-12, print_every_round=0)
theta_optim = numqi.optimize.minimize(model,**kwargs)
print(f'lower bound given by gradient descent = {theta_optim.fun}')
lower bound given by PPT = 0.3823670521387784
lower bound given by gradient descent = 0.3823671049768228
Both of them can detect the entanglement of the given state. However, we can go further by calculating the lower bound of $E_3(\rho)$ by gradient descent, which indicates the existence of higher entanglement in the given state.
dimA = 4
dimB = 4
tmp0 = [
[(0,0,1), (1,1,1), (2,2,1), (3,3,1)],
[(0,1,1), (1,2,1), (2,3,1), (3,0,1)],
[(0,2,1), (1,3,1), (2,0,1), (3,1,-1)],
]
matrix_subspace = np.stack([numqi.matrix_space.build_matrix_with_index_value(dimA, dimB, x) for x in tmp0])
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel((dimA, dimB), rank=2)
model.set_target(matrix_subspace)
kwargs = dict(theta0='uniform', num_repeat=3, tol=1e-12, print_every_round=0)
theta_optim = numqi.optimize.minimize(model,**kwargs)
print(f'lower bound of E_3 given by gradient descent = {theta_optim.fun}')
lower bound of E_3 given by gradient descent = 0.06557896267490304
Completely entangled subspace¶
Consider the following four three-qubit states
$$ \begin{aligned} \left|\psi_0\right\rangle&=|0\rangle_A|0\rangle_B|0\rangle_C, \\ \left|\psi_1\right\rangle&=|1\rangle_A|+\rangle_B|-\rangle_C, \\ \left|\psi_2\right\rangle&=|-\rangle_A|1\rangle_B|+\rangle_C, \\ \left|\psi_3\right\rangle&=|+\rangle_A|-\rangle_B|1\rangle_C. \end{aligned} $$ They form a UPB, which means no product state can be found in the orthogonal space of the subspace $\mathcal{S}$ spanned by them. In other words, $\mathcal{S}^{\bot}$ is a completely entangled subspace. We can use the gradient descent to calculate $E_2(\mathcal{S}^{\bot})$, which confirms the complete entanglement of the subspace.
rho_bes = numqi.entangle.load_upb('genshifts', 3, return_bes=True)[1]
EVL,EVC = np.linalg.eigh(rho_bes)
matrix_subspace = (EVC[:, EVL>1e-4]).reshape(2,2,2,-1).transpose(3,0,1,2)
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel((2,2,2), rank=1)
model.set_target(matrix_subspace)
kwargs = dict(theta0='uniform', num_repeat=3, tol=1e-12, print_every_round=0)
theta_optim = numqi.optimize.minimize(model,**kwargs)
print(f'gradient descent: {theta_optim.fun}')
gradient descent: 0.08144134645630785
The numerical result is really close to the analytical value $1-\frac{3\sqrt{6}}{8}$
Multipartite pure state¶
W state is a famous example of multipartite entangled state. It is defined as $$ \vert W \rangle = \frac{1}{\sqrt{3}}\left(\vert 001 \rangle + \vert 010 \rangle + \vert 100 \rangle\right). $$
We can calculate the value of $E_2(\vert W \rangle)$ and $E_3(\vert W \rangle)$ by gradient descent, which indicates the border rank of W state is 2.
num_qubit = 3
dim_list = [2]*num_qubit
Wstate = numqi.state.W(num_qubit).reshape(dim_list)
kwargs = dict(theta0='uniform', tol=1e-14, num_repeat=3, print_every_round=0, early_stop_threshold=1e-14)
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel(dim_list, rank=1)
model.set_target(Wstate)
theta_optim = numqi.optimize.minimize(model, **kwargs)
print(f'E_2 by gradient descent: {theta_optim.fun}')
model = numqi.matrix_space.DetectCanonicalPolyadicRankModel(dim_list, rank=2)
model.set_target(Wstate)
theta_optim = numqi.optimize.minimize(model, **kwargs)
print(f'E_3 by gradient descent: {theta_optim.fun}')
E_2 by gradient descent: 0.5555555555555554
E_3 by gradient descent: 8.530514072901951e-10