A Separability-Entanglement Classifier via Machine Learning

TL;DR

This paper introduces a machine learning-based classifier for quantum state entanglement, offering faster and more accurate identification compared to traditional methods, thus advancing quantum information processing.

Contribution

The paper presents a novel machine learning approach for quantum state separability classification, outperforming existing criteria in speed and accuracy.

Findings

Outperforms traditional methods in speed and accuracy

Effective for generic quantum states

Facilitates exploration of quantum entanglement

Abstract

The problem of determining whether a given quantum state is entangled lies at the heart of quantum information processing, which is known to be an NP-hard problem in general. Despite the proposed many methods such as the positive partial transpose (PPT) criterion and the k-symmetric extendibility criterion to tackle this problem in practice, none of them enables a general, effective solution to the problem even for small dimensions. Explicitly, separable states form a high-dimensional convex set, which exhibits a vastly complicated structure. In this work, we build a new separability-entanglement classifier underpinned by machine learning techniques. Our method outperforms the existing methods in generic cases in terms of both speed and accuracy, opening up the avenues to explore quantum entanglement via the machine learning approach.

Tables

Table 1. Table 1: Numerical results for approxmiating α tiles subscript 𝛼 tiles \alpha_{\text{tiles}} by α ​ ( 𝒞 , 𝒑 tiles ) 𝛼 𝒞 subscript 𝒑 tiles \alpha(\mathcal{C},\bm{p}_{\text{tiles}}) . Here, m 𝑚 m is the number of random extreme points for building 𝒞 𝒞 \mathcal{C} .

$m$	2000	5000	10000	20000	50000	100000
$\alpha (𝒞,{𝒑}_{\text{tiles}})$	0.5264	0.5868	0.6387	0.6759	0.7150	0.7459

Table 2. Table S2: The error rate of CHA 𝒞 m subscript 𝒞 𝑚 \mathcal{C}_{m} for two-qubit separable set and BCHA based on 𝒞 m subscript 𝒞 𝑚 \mathcal{C}_{m} for some critical m 𝑚 m .

$m$	1000	2000	3000	4000	5000
error of CHA (%)	8.55	6.01	4.85	4.05	3.60
error of BCHA (%)	3.03	1.97	1.47	1.17	1.15
$m$	6000	7000	8000	9000	10000
error of CHA (%)	3.25	2.95	2.76	2.64	2.55
error of BCHA (%)	1.01	0.75	0.79	0.71	0.65

Table 3. Table S3: Error rate of classifiers trained by different algorithms. The error rate is difficult to be reduced due to the lack of prior knowledge.

Method	Bagging	Boosting	SVM(rbf)	Decision Tree
Error (%)	12.03	14.8	8.4	23.3

Table 4. Table S4: The error rate of CHA 𝒞 m subscript 𝒞 𝑚 \mathcal{C}_{m} for two-qutrit separable set and BCHA based on 𝒞 m subscript 𝒞 𝑚 \mathcal{C}_{m} for some critical m 𝑚 m . Since we used 𝒞 10 5 subscript 𝒞 superscript 10 5 \mathcal{C}_{10^{5}} as the criterion, the error rate of 𝒞 10 5 subscript 𝒞 superscript 10 5 \mathcal{C}_{10^{5}} is 0 0 .

$m$	10000	20000	30000	40000	50000
error of CHA (%)	33.40	22.69	16.64	12.60	9.63
error of BCHA (%)	12.23	9.54	7.52	6.07	5.03
$m$	60000	70000	80000	90000	100000
error of CHA (%)	6.86	4.64	2.95	1.39	0
error of BCHA (%)	3.75	2.73	1.81	1.02	0

Equations

$$\mathcal{L}(h,{\mathcal{D}}_{\text{train}})=\frac{1}{|{\mathcal{D}}_{\text{train}}|}\sum _{({\mathbf{x}}_i,y_i)\in {\mathcal{D}}_{\text{train}}}\text{mathbbm}1(y_i\ne h({\mathbf{x}}_i)),$$

$$\max$$

$$\alpha \mathbf{p}=\sum _{i=1}^m{\lambda }_i{\mathbf{c}}_i,{\lambda }_i\ge 0,{\sum }_i{\lambda }_i=1.$$

$$\mathcal{L}(h,{\mathcal{D}}_{\text{train}})=\frac{1}{|{\mathcal{D}}_{\text{train}}|}\sum _{\left({\mathbf{x}}_i,{\alpha }_i,y_i\right)\in {\mathcal{D}}_{\text{train}}}\text{mathbbm}1(y_i\ne h({\mathbf{x}}_i,{\alpha }_i)).$$

$$s_{j,k}=E_{j,k}+E_{k,j}$$

$$a_{j,k}=-i\left(E_{j,k}-E_{k,j}\right)$$

$$d_l=\sqrt{\frac{2}{l(l+1)}}\left(\sum _{j=1}^l{E}_{j,j}-l{E}_{l+1,l+1}\right)$$

$$\text{tr}\left({\lambda }_i\right)=\text{tr}\left({\lambda }_i\mathbb{I}\right)=0$$

$$\text{tr}\left({\lambda }_i{\lambda }_j\right)=2{\delta }_{ij},$$

$$\rho =\frac{1}{n}\left(\mathbb{I}+\sqrt{\frac{n(n-1)}{2}}\mathbf{x}\cdot \vec{\lambda }\right),$$

$$x_i=\sqrt{\frac{n}{2(n-1)}}\text{tr}\left(\rho {\lambda }_i\right).$$

$${\Delta }_{\lambda }\left(d_1,\dots ,d_n\right)=C_{\lambda }\prod _{i=1}^n{d}_i^{-\lambda },$$

$${\mathcal{C}}_m:=\text{conv}\left(\left\{\text{0},\dots ,{\text{c}}_m\right\}\right)$$

$$\max$$

$$\alpha \text{p}=\sum _{i=0}^m{\lambda }_i{\text{c}}_i,$$

$${\lambda }_i\ge 0,{\sum }_i{\lambda }_i=1.$$

Full text

Born from pattern recognition, machine learning possesses the capability to make decisions without being explicitly programmed after learning from large amount of data. Beyond its extensive applications in industry, machine learning has also been employed to investigate physics-related problems in recent years. A number of promising applications have been proposed to date, such as the Hamiltonian learning Wiebe et al. (2014), automated quantum experiments generation Krenn et al. (2016), identification of phases and phase transition Schoenholz et al. (2016); van Nieuwenburg et al. (2017); Carrasquilla and Melko (2017), efficient representation of quantum many-body states Chen et al. (2017); Gao and Duan (2017), just to name a few. Nevertheless, there are yet a myriad of significant but hard problems in physics to be assessed, in which should machine learning provide more novel insights. For example, to determine whether a generic quantum state is entangled or not is a fundamental and NP-hard problem in quantum information processing Gurvits (2003), and machine learning is demonstrated to be exceptionally effective in tackling it as shown in this work.

As one of the key features in quantum mechanics, entanglement allows two or more parties to be correlated in a way that is much stronger than they can be in any classical way Horodecki et al. (2009). It also plays a key role in many quantum information processing tasks such as teleportation and quantum key distribution Nielsen and Chuang (2010). As a result, one question naturally arises: is there a universal criterion to tell if an arbitrary quantum state is separable or entangled? This is a typical classification problem, which remains of great challenge even for bipartite states. In fact, such an entanglement detection problem is proved to be NP-hard Gurvits (2003), implying that it is almost impossible to devise an efficient algorithm in complete generality.

Here, we focus on the task of detecting bipartite entanglement. Consider a bipartite system $AB$ with the Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, where $\mathcal{H}_{A}$ has dimension $d_{A}$ and $\mathcal{H}_{B}$ has dimension $d_{B}$, respectively. A state $\rho_{AB}$ is separable if it can be written as a convex combination $\rho_{AB}=\sum_{i}\lambda_{i}\rho_{A,i}\otimes\rho_{B,i}$ with a probability distribution $\lambda_{i}\geq 0$ and $\sum_{i}\lambda_{i}=1$. Here $\rho_{A,i}$ and $\rho_{B,i}$ are density operators acted on $\mathcal{H}_{A}$, $\mathcal{H}_{B}$ respectively. Otherwise, $\rho_{AB}$ is entangled Werner (1989); Gühne and Tóth (2009). To date, many criteria have been proposed to detect bipartite entanglement, each with its own pros and cons. For instance, the most famous criterion is the positive partial transpose (PPT) criterion, saying that a separable state must have PPT; however, it is only necessary and sufficient when $d_{A}d_{B}\leq 6$ Peres (1996); Horodecki et al. (1996). Another widely used one is the $k$-symmetric extension hierarchy Navascués et al. (2009); Brandão et al. (2011), which is presently one of the most powerful criteria, but hard to compute in practice due to its exponentially growing complexity with $k$ S .

In this work, we employ the machine learning techniques to tackle the bipartite entanglement detection problem by recasting it as a learning task, namely we attempt to construct a separability-entanglement classifier. Due to its renowned effectiveness in pattern recognition for high-dimensional objects, machine learning is a powerful tool to solve the above problem. In particular, a reliable separability-entanglement classifier in terms of speed and accuracy is constructed via the supervised learning approach. The idea is to feed our classifier by a large amount of sampled trial states as well as their corresponding class labels (separable or entangled), and then train the classifier to predict the class labels of new states that it has not encountered before. It is worthy stressing that, there is also a remarkable improvement with respect to universality in our classifier compared to the conventional methods. Previous methods only detect a limited part of the state space, e.g. different entangled states often require different entanglement witnesses. In contrast, our classifier can handle a variety of input states once properly trained, as shown in Fig. 1.

Supervised learning – The bipartite entanglement detection problem can be formulated as a supervised binary classification task. Following the standard procedure of supervised learning Mohri et al. (2012); Shalev-Shwartz and Ben-David (2014), the feature vector representation of the input objects (states) in a bipartite system AB is first created. Indeed, any quantum state $\rho$, as a density operator acting on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ can be represented as a real vector $\bm{x}$ in $\mathscr{X}=\mathbb{R}^{d_{A}^{2}d_{B}^{2}-1}$, which is due to the fact that $\rho$ is Hermitian and of trace $1$ (see Supplementary Material S ). In the machine learning language, we refer $\bm{x}$ as the feature vector of $\rho$ and $\mathscr{X}$ the feature space.

Next, a dataset of training examples is produced, with the form $\mathcal{D}_{\text{train}}=\left\{(\bm{x}_{1},y_{1}),...,(\bm{x}_{n},y_{n})\right\}$, where $n$ is the size of the set, $\bm{x}_{i}\in\mathscr{X}$ is the $i$-th sample, and $y_{i}$ is its corresponding label signifying which class it belongs to, i.e., $y_{i}$ equals to 1 if $\bm{x}_{i}$ is entangled or $-1$ otherwise. When $d_{A}d_{B}\leq 6$, the labeling process can be directly computed via the PPT criterion. For higher-dimensional cases, we attempt to estimate the labels by convex hull approximation, which we will describe later. The task is to analyze these training data and produce an inferred classifier that predicts the unknown class labels for generic new input states.

Explicitly, the aim of supervised learning is to infer a function (classifier) $h:\mathscr{X}\to\left\{-1,1\right\}$ among a fixed class of functions $\mathscr{H}$ such that $h$ is expected to be close to the true decision function. One basic approach to choose $h$ is the so-called empirical risk minimization, which seeks the function that best fits the training data among the class $\mathscr{H}$. In particular, to evaluate how well $h$ fits the training data $\mathcal{D}_{\text{train}}$, a loss function is defined as

[TABLE]

where $\mathbbm{1}(\cdot)$ is the truth function of its arguments. For a generic new input test dataset $\mathcal{D}_{\text{test}}$ that contains previously unseen data, function $\mathcal{L}(h,\mathcal{D}_{\text{test}})$ gives a quantification of the generalization error from $\mathcal{D}_{\text{train}}$ to $\mathcal{D}_{\text{test}}$.

Numerous supervised learning algorithms have been developed, each with its strength and weakness. These algorithms, which have distinct choices of class $\mathscr{H}$, include support vector machine (SVM) Cortes and Vapnik (1995), decision trees Breiman et al. (1984), bootstrap aggregating Zhou (2012), and boosting Schapire (2003), etc. We have applied these algorithms to the separability problem directly, but neither of them provided an acceptable accuracy, which is mainly due to the lack of prior knowledges for training, e.g., the geometric shape of the set of separable states $\mathcal{S}$. Taking the kernel SVM approach Cortes and Vapnik (1995) as an example, it uses a kernel function to map data from the original feature space to another Hilbert space, and then finds a hyperplane in the new space to split the data into two subclasses. It turns out that using common kernels such as radial basis function and polynomials, the error rate on the test dataset is always around $10\%$ (see Supplementary Material S for details). This suggests that the boundary of $\mathcal{S}$ is too complicated to be portrayed by manifolds with ordinary shapes.

Convex hull approximation – The above discussions suggest that it is desirable to examine the detailed geometric shape of $\mathcal{S}$ in advance. One well-known approach is to approximate $\mathcal{S}$ from outsize via $k$-symmetric extendible set $\Theta_{k}$, where $\Theta_{k}\supset\Theta_{k+1}$ and $\Theta_{k}$ converges exactly to $\mathcal{S}$ as $k$ goes to infinity Doherty et al. (2002). Unfortunately, it is impractical to compute the boundary of $\Theta_{k}$ for large $k$, while it is still far from approximating $\mathcal{S}$ for small $k$ S .

However, it is much easier to approximate $\mathcal{S}$ from inside, since $\mathcal{S}$ is a closed convex set, and its extreme points are exactly all the separable pure states, which can be straightforwardly parameterized and generated numerically. We randomly sample $m$ separable pure states $\bm{c}_{1},...,\bm{c}_{m}\in\mathscr{X}$ to form a convex hull $\mathcal{C}\coloneqq\operatorname{conv}(\left\{\bm{c}_{1},...,\bm{c}_{m}\right\})$. $\mathcal{C}$ is said to be a convex hull approximation (CHA) of $\mathcal{S}$, with which we can approximately tell whether a state $\rho$ is separable or not by testing if its feature vector $\bm{p}$ is in $\mathcal{C}$. This is equivalent to determining whether $\bm{p}$ can be written as a convex combination of $\bm{c}_{i}$ by solving the following linear programming:

[TABLE]

Here $\alpha=\alpha(\mathcal{C},\bm{p})$ is a function of $\mathcal{C}$ and $\bm{p}$. If $\alpha(\mathcal{C},\bm{p})\geq 1$, $\bm{p}$ is in $\mathcal{C}$ and thus $\rho$ is separable; otherwise, $\rho$ is highly possible to be an entangled state. In principle, $\mathcal{C}$ will be a more accurate CHA of $\mathcal{S}$ if we construct $\mathcal{C}$ with more extreme points. We test the error rate of CHA on a set of $2\times 10^{4}$ random two-qubit states, which is sampled under a specified distribution  S and labeled by PPT criterion. The results are shown by the blue curve in Fig. 3(c), where the error rate decreases quickly to $3\%$ when the number of extreme points $m$ increases to $10^{4}$.

However, we can not directly test the accuracy of CHA on generic two-qutrit states, since PPT criterion is no longer sufficient for detecting separability. To illustrate the power of CHA beyond the PPT criterion, we use a specific example that is previously well-studied. Consider a set of two-qutrits pure states $\{|v_{1}\rangle,\dots,|v_{5}\rangle\}$ that form the well-known unextendible product basis Bengtsson and Zyczkowski (2007), where $|v_{1}\rangle=(|00\rangle-|01\rangle)/\sqrt{2}$, $|v_{2}\rangle=(|21\rangle-|22\rangle)/\sqrt{2}$, $|v_{3}\rangle=(|02\rangle-|12\rangle)/\sqrt{2}$, $|v_{4}\rangle=(|10\rangle-|20\rangle)/\sqrt{2}$, and $|v_{5}\rangle=(|0\rangle+|1\rangle+|2\rangle)^{\otimes 2}/3$. It is known that $\rho_{\text{tiles}}=(\mathbb{I}-\sum_{i=1}^{5}|v_{i}\rangle\langle v_{i}|)/4$ is an entangled state with PPT Bennett et al. (1999). Due to the fact that $\mathcal{S}$ is convex and closed, there must exist a unique critical point $\alpha_{\text{tiles}}\in[0,1)$ such that $\alpha_{\text{tiles}}\rho+(1-\alpha_{\text{tiles}})\mathbb{I}/(d_{A}d_{B})$, the probabilistic mixture of $\rho$ and the maximally-mixed state $\mathbb{I}/(d_{A}d_{B})$, is on the boundary of $\mathcal{S}$. Ref. Johnston (2014) compared the effectiveness of various separability criteria, and concluded that $\alpha_{\text{tiles}}\in(0.5643,0.8649]$. Note that for a CHA $\mathcal{C}$, $\alpha(\mathcal{C},\bm{p}_{\text{tiles}})$ actually provides a lower bound approximation of $\alpha_{\text{tiles}}$, where $\bm{p}_{\text{tiles}}$ is the feature vector of $\rho_{\text{tiles}}$. Now we apply CHA and attempt to improve the lower bound of $\alpha_{\text{tiles}}$, and the result is shown in Table 1.

We find that the lower bound of $\alpha_{\text{tiles}}$ has been raised to $0.7459$. However, in Table 1, the value of $\alpha(\mathcal{C},\bm{p}_{\text{tiles}})$ has not converged yet. To reach the convergence of $\alpha(\mathcal{C},\bm{p}_{\text{tiles}})$, we have to enlarge $\mathcal{C}$ by adding more extreme points. However, note that the point $\alpha(\mathcal{C},\bm{p}_{\text{tiles}})$ lies on a part of the boundary of $\mathcal{C}$, which is the intersection of a hyperplane and $\mathcal{C}$. Let $\bm{c}_{i_{1}},\ldots,\bm{c}_{i_{D}}$ be the extreme points of $\mathcal{C}$ that lie on the hyperplane as well. Clearly, if we enlarge $\mathcal{C}$ by sampling the separable pure states that are near $\bm{c}_{i_{1}},\ldots,\bm{c}_{i_{D}}$ rather than sampling uniformly over the whole set of separable pure states, it will boost the value of $\alpha(\mathcal{C},\bm{p}_{\text{tiles}})$ more effectively.

Subsequently, we refine CHA as an iterative algorithm  S , with the idea shown in Fig. 2(a). The iterative algorithm gives the result $\alpha_{\text{tiles}}>0.8648$. As the upper bound of $\alpha_{\text{tiles}}$ is $0.8649$  Johnston (2014), we can explicitly conclude that $\alpha_{\text{tiles}}\approx 0.8649$. It is worthy emphasizing that the algorithm also gives the critical point for a generic entangled state with small error, and detects the separability for generic separable states S .

Combining CHA and supervised learning – There is yet a noticeable drawback of the above CHA approach from the perspective of the tradeoff between the accuracy and time consumption. Boosting the accuracy means adding additional extreme points to enlarge the convex hull, which leads to more time costs to determine if a point is inside the enlarged convex hull or not. To overcome this, we combined CHA with supervised learning, as machine learning has the power to speed up such computations.

To design a learning process that is suitable for our problem, for each state $\rho$ with feature vector $\bm{p}$, we extend the feature vector as $(\bm{p},\alpha(\mathcal{C},\bm{p}))$ in order to encode the geometric information of the CHA $\mathcal{C}$ into the dataset. In this manner, the training dataset is written as $\mathcal{D}_{\text{train}}=\left\{(\bm{x}_{1},\alpha_{1},y_{1}),...,(\bm{x}_{n},\alpha_{n},y_{n})\right\}$, where $\alpha_{i}=\alpha(\mathcal{C},\bm{x}_{i})$. A classifier $h$ is now a binary function defined on $\mathscr{X}\times\mathbb{R}$, and the loss function of a classifier $h$ is then redefined as

[TABLE]

Subsequently, we employ a standard ensemble learning approach Zhou (2012) to train a classifier with training data $\mathcal{D}_{\text{train}}$, described as follows.

The essential idea of ensemble learning is to improve the predictive performance by combining multiple classifiers into a committee. Even though the prediction from each constituent might be poor, the combined classfier could often still perform excellent. For the binary classification problem, we can train different classifiers to give their respective binary votes for each prediction, and use the majority rule to choose the value which receives more than half votes as the final answer, see Fig. 2(c) for a schematic diagram.

Here, we choose bootstrap aggregating (bagging) Breiman (1996) as our training ensemble algorithm. In each run, a training subset is randomly drawn from the whole set $\mathcal{D}_{\text{train}}$, and a model is trained from the training subset using another learning algorithm, e.g., decision trees learning. We repeat the process for $L=100$ times and obtain $L$ different models, which are finally combined together as the committee. Since $\alpha_{i}=\alpha(\mathcal{C},\bm{x}_{i})$ contains the geometric information of CHA $\mathcal{C}$, our method is indeed a combination of bagging and CHA. We call this combined method BCHA.

The computational cost contains two parts: the cost of computing $\alpha_{i}$ via linear programming, and the time of computing each constituent in the committee. The latter cost is much smaller than the former. Therefore, by using a convex hull of much smaller size and implementing a bagging algorithm, a significant boost in terms of accuracy is anticipated if the total computational cost is fixed. For the two-qubit case, we have demonstrated such a remarkable boost of accuracy in our BCHA classifier, as shown in Fig. 3(c), where the advantages of the BCHA classifier in terms of both accuracy and speed are shown.

We further extend the classifier to the two-qutrit scenario. Unlike the two-qubit case, the critical question now is how to set an appropriate criterion to evaluate whether the classifier is working correctly, since PPT criterion is not sufficient for detecting separability in two-qutrit systems. As the convex hull is capable of approximating the set of separable states $\mathcal{S}$ to an arbitrary precision, we use $10^{5}$ random separable pure states as extreme points to form the hull, and assumed it to be the true $\mathcal{S}$. The learning procedure is analogous to the one used for two qubits. Figure 3(d) shows the accuracy of the BCHA classifier compared to that of the sole CHA approach, Similar as the two-qubit case, the BCHA classifier shows clear advantage in terms of both accuracy and speed in comparison with the sole CHA method.

Conclusion – In summary, we study the entanglement detection problem via the machine learning approach, and build a reliable separability-entanglement classifier by combining supervised learning and the CHA method. Compared to the conventional criteria for entanglement detection, our method can classify an unknown state into the separable or entangled category more precisely and rapidly. The classifier can be extended to higher dimensions in principle, and the developed techniques in this work would also be incorporated in future entanglement-engineering experiments. We anticipate that our work would provide new insights to employ the machine learning techniques to deal with more quantum information processing tasks in the near future.

Appendix A A. Generalized Gell-Mann Matrices

To represent a $n$-by-$n$ density matrix $\rho$ as a real vector $\bm{x}$ in $\mathbb{R}^{n^{2}-1}$, we can find a Hermitian orthogonal basis that contains identity such that $\rho$ can be expanded in such a basis with real coefficients. For example, the Pauli basis is a commonly used one. In our numerical tests, we take the generalized Gell-Mann matrices and the identity as the Hermitian orthogonal basis. In this section, we recall the definition of the generalized Gell-Mann matrices, which is shown in Bertlmann and Krammer (2008).

Let $\left\{|1\rangle,\ldots,|n\rangle\right\}$ be the computatioal basis of the $n$-dimensional Hilbert space, and $E_{j,k}=|j\rangle\langle k|$. We now define three collections of matrices. The first collection is symmetric:

[TABLE]

for $1\leq j<k\leq n$. The second collection is antisymmetric:

[TABLE]

for $1\leq j<k\leq n$. The last collection is diagonal:

[TABLE]

for $1\leq l\leq n-1$.

The generalized Gell-Mann matrices are elements in the set $\left\{\lambda_{i}\right\}=\left\{s_{j,k}\right\}\cup\left\{a_{j,k}\right\}\cup\left\{d_{l}\right\}$, which gives a total of $n^{2}-1$ matrices. We can easily check that

[TABLE]

and

[TABLE]

which implies that $\left\{\lambda_{i}\right\}\cup\left\{\mathbb{I}\right\}$ forms an orthogonal basis of observables in $n$-dimensional Hilbert space.

For every $n$-by-$n$ density matrix $\rho$, $\rho$ can be expressed as a linear combination of $\lambda_{i}$ and $\mathbb{I}$ as follows:

[TABLE]

where $\bm{x}=(x_{1},x_{2},\dots,x_{n^{2}-1})\in\mathbb{R}^{n^{2}-1}$ satisfies

[TABLE]

Appendix B B. The set of $k$-extendible states

In this section, we recall facts regarding $k$-extendible states and its relationship to separability.

A bipartite state $\rho_{AB}$ is said to be $k$-symmetric extendible if there exists a global state $\rho_{AB_{1}\ldots B_{k}}$ whose reduced density matrices $\rho_{AB_{i}}$ are equal to $\rho_{AB}$ for $i=1,\ldots,k$. The set of all $k$-extendible states, denoted by $\Theta_{k}$, is convex with a hierarchy structure $\Theta_{k}\supset\Theta_{k+1}$. Moreover, when $k\rightarrow\infty$, $\Theta_{k}$ converges exactly to the set of separable states Doherty et al. (2002).

The $\Theta_{k}$ is known to be closely related to the ground state of some $(k+1)$-body Hamiltonians Zeng et al. (2015). To be more precise, consider a $2$-local Hamiltonian $H$ of a $(k+1)$-body system with Hilbert space $\mathbb{C}^{d_{A}}\bigotimes_{i=1}^{k}\mathbb{C}^{d_{B_{i}}}$ of dimension $d_{A}d_{B}^{k}$, as given in the following form $H=\sum_{i=1}^{k}H_{AB_{i}}.$ Here $H_{AB_{i}}$ is any Hermitian operator acting nontrivially on particles $A$ and $B_{i}$, and trivially on other $k-1$ parties. In other words, we will have $H_{AB_{1}}=H_{AB}\otimes I_{2,\ldots,k}$ ($I_{2,\ldots,k}$ is the identity operator of $B_{2},\ldots,B_{k}$), and given the symmetry of $B_{i}$s, we can always write the nontrivial action of $H_{AB_{i}}$ on $\mathbb{C}^{d_{A}}\bigotimes_{i=1}^{k}\mathbb{C}^{d_{B_{i}}}$ in terms of some $H_{AB}$ acted on $d_{A}d_{B}$-dimensional Hilbert space.

For any given $H$, denote its normalized ground state by $|\psi_{g}\rangle\in\mathbb{C}^{d_{A}}\bigotimes_{i=1}^{k}\mathbb{C}^{d_{B_{i}}}$, and $\rho_{g}=|\psi_{g}\rangle\langle\psi_{g}|$. Then the extreme points of $\Theta_{k}$ are given by the marginals of $\rho_{g}$ on particles $AB_{i}$, which are the same for any $i$. Denote this marginal by $\rho_{H}$, since it is completely determined by $H$.

To generate random extreme points of $\Theta_{k}$, we will need to first parametrize them. Denote $\{O_{AB}^{lm}\}$ as set of orthonormal Hermitian basis for operators on $H_{A}\otimes H_{B_{j}}$ (see section A), then we can always write $H_{AB}=\sum_{lm}a_{lm}O_{AB}^{lm},$ with parameters $a_{lm}$. Without loss of generality, we assume $O_{AB}^{00}=I$, and we will assume $a_{00}=0$, so there are only $d_{A}^{2}d_{B}^{2}-1$ terms in the sum. Since $H_{AB}$ is a Hermitian matrix, $a_{lm}$ can be chosen as real, and we can further require that $\sum_{l,m}a_{lm}^{2}=1.$ Consequently, $\rho_{H}$ will be a point in $\mathbb{R}^{d_{A}^{2}d_{B}^{2}-1}$, which is parametrized by $\{a_{lm}\}$. And the coordinate of $\rho_{H}$ are explicitly given by $b_{lm}=\operatorname{tr}(\rho_{H}O_{AB}^{lm})$.

Also, each $H_{AB}$ gives an entanglement witness. The ground state energy of $H_{AB}$ is given by $E_{0}=\langle\psi_{g}|H|\psi_{g}\rangle=\sum_{lm}a_{lm}b_{lm}/k.$ For any density matrix $\rho_{AB}$, if $\operatorname{tr}(\rho_{AB}H_{AB})<E_{0}$, then $\rho_{AB}$ has no $k$-symmetric extension, hence is surely entangled.

Since the dimension of $H$ grows exponentially with $k$, to generate these extreme points for $\Theta_{k}$ becomes hard when $k$ increases. In practice, we can generate the extreme points of $\Theta_{k}$ for $k=12$ and $d_{A}=d_{B}=2$. However, as depicted in Fig. S4, there is still a large gap between the separable boundary and the $k$-extension boundary $k$.

More general properties on $k$-extendability and its relationship to the quantum marginal problem can be found in Refs.  Klyachko (2006); Coleman (1963); Erdahl (1972); Liu (2006); Liu et al. (2007); Wei et al. (2010).

Appendix C C. Generating Random Density Matrices

Since our aim is to determine the separability of generic bipartite states, we require a bunch of random density matrices with full rank to test the performance of our approaches. In our numerical tests, we sample random density matrices under the probability measure $\mu=\nu\times\Delta_{\lambda}$, where $\nu$ is the uniform distribution on $U(n)$ according to Haar measure, $\Delta_{\lambda}$ is the Dirichlet distribution on the simplex $\sum_{i=1}^{n}d_{i}=1$. The probability density function of Dirichlet distribution is

[TABLE]

where $\lambda>0$ is a parameter and $C_{\lambda}$ is the normalization constant. Since every density matrix is unitarily similar to a real diagonal density matrix, $\mu$ is a probability measure on the set of all density matrices. Such a probability measure is discussed in Ref. Życzkowski (1999), section II.A. We implemented the sampling on $\nu$ via directly calling the function RandomUnitary in Johnston (2016). The entire implementation is in the code RandomState.m on the website of QMLab qml (2017).

For the 2-qubit case, we set $\lambda=1/2$ and generate $5\times 10^{4}$ random quantum states, which are put in the file 2x2rdm.mat. We have found that $35\%$ of the states are PPT states, i.e., separable states, which is consistent with the result shown in Życzkowski (1999).

For the 2-qutrit case, we also set $\lambda=1/2$. As shown in Życzkowski (1999), only $2.2\%$ of the random states are PPT states when $\lambda=1/2$. However, our main interest is determining whether a PPT state is entangled. Thus, we reject all the states with negative partial transpose while sampling, and obtain $2\times 10^{4}$ PPT states eventually. These states are put in the file 3x3rdm.mat on qml (2017). We can verify that at least $66.24\%$ of the PPT states are separable under the probability measure we have chosen, using the convex hull approximation, which will be discussed in section E.

We determine whether a state is a PPT state via IsPPT in Johnston (2016).

Appendix D D. Testing CHA and BCHA

To approximate the set of separable states $\mathcal{S}$ with a convex hull $\mathcal{C}$, we generate a bunch of extreme points of $\mathcal{S}$, i.e., random separable pure states in $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ in a straightforward way. The procedure for each time of sampling is demonstrated as follows:

Sample a state vector $|\psi_{A}\rangle\in\mathcal{H}_{A}\cong\mathbb{C}^{d_{A}}$ from uniform distribution on the unit hypersphere in $\mathbb{C}^{d_{A}}$, according to Haar measure  Johnston (2016). 2. 2.

Sample another state vector $|\psi_{B}\rangle\in\mathcal{H}_{B}\cong\mathbb{C}^{d_{B}}$ from uniform distribution on the unit hypersphere in $\mathbb{C}^{d_{B}}$, according to Haar measure. 3. 3.

Return $|\psi_{A}\rangle|\psi_{B}\rangle$.

We execute the above procedure for $M$ times to gain $M$ extreme points $\textbf{c}_{1},\ldots,\textbf{c}_{M}$. Let

[TABLE]

for $m=1,\ldots,M$. It is easy to see that $\mathcal{C}_{m}\subseteq\mathcal{C}_{m+1}$ for $m=1,\ldots,M-1$. Recall that we can decide whether a point p is in $\mathcal{C}_{m}$ by solving the following linear programming

[TABLE]

If $\alpha\geq 1$, p is in $\mathcal{C}_{m}$ and thus separable; otherwise, it is possibly an entangled state. The solver for the linear programming 4 is implemented in CompAlpha.m on qml (2017).

For the two-qubit case, we sample $M=10^{4}$ extreme points, which is saved in the file 2x2extreme.mat on qml (2017). We split the data in 2x2rdm.mat into two, one for training BCHA and the other for testing both CHA and BCHA. To compare the performance of CHA and BCHA, we test the error rate of the CHA $\mathcal{C}_{m}$ and BCHA based on $\mathcal{C}_{m}$ on the test dataset. The result is shown in table S2.

We also apply different supervised learning algorithms with the same training and test dataset, without combining CHA. The result is shown in Table S3.

For the two-qutrit case, we sample $M=10^{5}$ extreme points, which is saved in the file 3x3extreme.mat on qml (2017). It can be verified that $66.24\%$ of the PPT random states in 3x3rdm.mat are in the convex hull $\mathcal{C}_{10^{5}}$, which implies that at least $66.24\%$ of the PPT random states are separable. We used $\mathcal{C}_{10^{5}}$ as the criterion for separability, i.e., regard $\mathcal{C}_{10^{5}}$ as the true separable set. Similar to the two-qubit case, we also tested the accuracy of CHA $\mathcal{C}_{m}$ as well as BCHA based on $\mathcal{C}_{m}$. The result is shown in table S4.

Appendix E E. Iterative Algorithm for Computing the Critical Point

Recall that for an entangled state $\rho$, there exists a critical point $\alpha_{\rho}$ such that $\alpha\rho+(1-\alpha)\mathbb{I}/(d_{A}d_{B})\ \;(0\leq\alpha\leq 1)$ is separable when $\alpha\leq\alpha_{\rho}$ and entangled when $\alpha>\alpha_{\rho}$. Based on CHA, we developed an iterative algorithm for approximating $\alpha_{\rho}$ in a more efficient way, which is shown as follows:

Randomly sample $1000$ extreme points and form a convex hull $\mathcal{C}$. Let p be the feature vector of $\rho$. Set $\epsilon=1$, $\gamma=0.95$. 2. 2.

Update $\alpha_{\rho}\leftarrow\alpha(\mathcal{C},\textbf{p})$. 3. 3.

Suppose now $\mathcal{C}:=\text{conv}\left(\left\{\textbf{c}_{1},\ldots,\textbf{c}_{m}\right\}\right)$, and $\alpha_{\rho}\textbf{p}=\sum_{i}\lambda_{i}\textbf{c}_{i}$. Let $\textbf{c}_{i_{1}},\ldots,\textbf{c}_{i_{D}}$ be the extreme points such that $\lambda_{i_{k}}>0$. Set $\mathcal{C}\leftarrow\text{conv}\left(\left\{\textbf{c}_{i_{1}},\ldots,\textbf{c}_{i_{D}}\right\}\right)$. 4. 4.

For each $k=1,\ldots,D$, suppose $\textbf{c}_{i_{k}}$ is the feature vector of $|a_{k}\rangle|b_{k}\rangle$. We randomly generate two Hermitian operators $H_{1}\in\text{End}(\mathcal{H}_{A})$, $H_{2}\in\text{End}(\mathcal{H}_{B})$ such that $\left\lVert H_{1}\right\rVert_{2}=1$ and $\left\lVert H_{2}\right\rVert_{2}=1$. Let $\xi$ be a random number in $[0,\epsilon]$. Set $|a_{k}^{\prime}\rangle|b_{k}^{\prime}\rangle=\left(e^{i\xi H_{1}}\otimes e^{i\xi H_{2}}\right)|a_{k}\rangle|b_{k}\rangle$. Set $\mathcal{C}\leftarrow\text{conv}\left(\left\{\mathcal{C},\textbf{c}_{k}^{\prime}\right\}\right)$, where $\textbf{c}_{k}^{\prime}$ as the feature vector of $|a_{k}^{\prime}\rangle|b_{k}^{\prime}\rangle$. 5. 5.

Set $\epsilon\leftarrow\gamma\epsilon$ and back to step 2.

What step 4 does is sampling in the neighborhood of $\textbf{c}_{i_{k}}$. In practice, we repeat step 4 for $10$ times to get a bunch of neighbors. The detailed implementation is in the code CriticalPoint.m on  qml (2017).