# How to find a separable decomposition of a separable density matrix?

If I have a bipartite system of two qubits $$A$$ and $$B$$, and the density matrix $$\rho$$ is separable, how do I decompose it into its separable parts?

That is, give $$\rho$$, expand it as follows:

$$\rho = \sum_{i=0}^Np_i \ \rho_{i}^A \otimes \rho_{i}^B$$

Where $$0 \le p_i \le 1$$ and $$\rho^{A,\ B}_i$$ are density matrices on the two subsystems $$A$$ and $$B$$ for some $$N$$.

Cross-posted on qc.SE

• This is a hard problem. (NP-hard -- in the dimension of the spaces involved, IIRC.) Nov 6, 2022 at 22:52
• related: physics.stackexchange.com/q/399675/58382 (and posts on qc.SE linked in the qc.SE version of this question)
– glS
Nov 16, 2022 at 14:13

## 1 Answer

I also found myself interested in this problem and what seemed like a lack of discussion on how to tackle it (as pointed out in the comments it is, computationally, a very hard problem in general).

In the literature a most significant contribution came from Lewenstein and Sanpera in 1998 with a Letter proving that for any density matrix $$\rho$$ and set $$V$$ of product vectors belonging to the range (i.e. rowspace) of $$\rho$$ there exists a best separable approximation (BSA) of the matrix in the sense $$\rho=\rho_s^*+ \delta \rho$$ where $$\rho_s^*$$ is separable and the trace $$\text{Tr}\delta \rho$$ of the entangled part is minimised.

They also provided an explicit construction for the separable part. The idea is to generate a large set $$V$$ of projectors $$\lvert \alpha \rangle \langle \alpha \rvert$$ randomly chosen from $$\text{Range}(\rho)$$ and seek a decomposition $$\rho_s^*=\sum_{\alpha\in V} \Lambda_\alpha \lvert \alpha \rangle \langle \alpha \rvert$$ where the $$\Lambda_\alpha \geq 0$$ are pairwise maximised: $$(\Lambda_{\alpha_1} \lvert \alpha_1 \rangle \langle \alpha_1 \rvert + \Lambda_{\alpha_2} \lvert \alpha_2 \rangle \langle \alpha_2 \rvert)$$ should be the largest contribution that can subtracted from $$\rho$$ whilst retaining the non-negativity of the difference $$(\dagger)$$.

Lewenstein and Sanpera include explicit expressions for the $$\Lambda_\alpha$$, but in practice resort to trial-and-error (i.e. brute force) to sidestep potential issues of convergence here. Below is the one figure from their Letter showing the trace of the residual (non-separable) part when applying this approach to a $$2\otimes2$$ Werner state $$\rho_w(x)$$ known to be separable for values of the parameter $$x\leq 1/3$$. The upshot is that, provided the number of product vectors chosen for $$V$$ is sufficiently large, the result does not depend on this number and the true BSA is obtained (the non-zero part seen for $$x\leq 1/3$$ is a result of the numerical precision set to $$10^{-4}$$ here).

That's the general problem; following a brief flurry of activity after the Lewenstein-Sanpera result (the "Lewenstein-Sanpera" decomposition) I could not find much further work for arbitrary density matrices. However, there has been a great deal of progress for specific cases (Edit: as kindly pointed out in the comments, this includes extensive work in the domain semidefinite programming for bipartite states and hardness). In particular the problem you describe, that for two qubits, does in fact have a closed form solution.

The clearest presentation of this I could find is by Jafarizadeh et al. (2005) who rather elegantly show that the problem can be formulated as one of semidefinite programming and proceed to derive the best separable approximation of a generic two qubit density matrix using the Wootters basis $$\{\lvert x_i \rangle\}$$ for $$\rho$$ $$\displaystyle({\dagger\dagger})$$. In the case $$\rho$$ is separable this reduces to $$\rho=\rho_s$$ where \begin{align*} \rho_s = ( \lambda_2+\lambda_3+\lambda_4 ) \lvert x_1' \rangle \langle x_1' \rvert + \lambda_2 \lvert x_2' \rangle \langle x_2' \rvert + \lambda_3 \lvert x_3' \rangle \langle x_3' \rvert + \lambda_4 \lvert x_4' \rangle \langle x_4' \rvert \end{align*} is in fact a separable matrix. Here $$\lambda_i$$ are the squares roots of the eigenvalues in decreasing order of the non-Hermitian matrix $$\rho \tilde{\rho}$$, with $$\tilde{\rho} = (\sigma_y \otimes \sigma_y) \rho^* (\sigma_y \otimes \sigma_y),$$ and $$\lvert x_i'\rangle=\lvert x_i \rangle /\sqrt{\lambda_i}$$ are the Wootters states scaled by the square root of these values.

$$(\dagger)$$ e.g. For a single $$\Lambda_\alpha$$ this means $$\rho- \Lambda_\alpha \lvert \alpha \rangle \langle \alpha \rvert \geq 0$$ and for every $$\epsilon\geq0$$, $$\rho-\Lambda_\alpha \lvert \alpha \rangle \langle \alpha \rvert$$ is not positive definite.

$$(\dagger\dagger)$$ This is a (non-separable) decomposition $$\rho=\sum_i \lvert x_i \rangle \langle x_i \rvert$$ shown by Wootters (also in 1998) to always exist.

• Your answer sounds like no progress has been made on this in 20 years! Feb 11, 2023 at 18:37
• Quite surprising isn't it! I couldn't find anything on the general problem citing the Lewenstein-Sanpera letter in the last 15 years. I am not familiar with the topic but it seems like more robust and efficient algorithms for the BSA have been developed using convex programming (see doi.org/10.1016/j.laa.2006.08.026 using a Frank-Wolfe style algorithm). It would be interesting for an expert to weigh in on this.
– pip
Feb 12, 2023 at 6:58
• My comment was more intended as a criticism -- there has definitely been work on this after the late 90ies (though I am not an expert, so it would take some time to find it). Just two quick pointers -- there has been work on hardness, e.g. arxiv.org/abs/0810.4507 (and the work of Gurvits quoted there), and a substantial body of work on SDP relaxations and other algorithms for the separability problem, e.g.link.springer.com/article/10.1007/s00220-019-03382-y, arxiv.org/abs/1504.01029, arxiv.org/abs/1011.2751 (this is really just the result ... Feb 12, 2023 at 12:33
• ... of some quick searches, and by no means representative). Feb 12, 2023 at 12:33
• That's great, thanks! (There is something to be said about the fastest way to get an answer for your question on the internet is to first answer it incorrectly...). Are the results you link for bipartite constructions only? Shortly after posting I did come across doi.org/10.1016/j.laa.2006.08.026, but did not feel in a position to comment on it (I have no education in convex optimisation). It would be great if someone qualified could summarise these works for us. In the meantime, I have added an (Edit:) to hopefully point out the limitations of my answer.
– pip
Feb 12, 2023 at 17:42