How to factorize a multi-qubit quantum state?

Question

Let's say I have an arbitrary n-qubit quantum state $|\psi\rangle$. Is there a simple way to compute the states of individual qubits (or sets of qubits if they are entangled) that compose $|\psi\rangle$? E.g. if $n=3$ how can I factorize $|\psi\rangle$ into $|\psi\rangle = |q_0\rangle\otimes |q_1\rangle \otimes |q_2\rangle$ or if let's say $|q_0\rangle$ and $|q_1\rangle$ are entangled $|\psi\rangle=|q_0q_1\rangle\otimes|q_2\rangle$? Or as the case may be, how can I verify it's not factorizable? It's an easy task for the case of 2 qubits, but in general case it looks like I need to test all possible combinations of entangled subsets?

Answer 1

What you ask is even harder than asking whether a state is separable. The latter is considered a complex problem, sometimes NP-hard. See, e.g., https://en.wikipedia.org/wiki/Separable_state for a quick overview. So, to answer your question, "Is there a simple way...", probably no, given that you assume an arbitrary state $|\psi\rangle$.

Answer 2

You ask two slightly different (but somewhat related) questions about computing "the states of individual qubits" vs. finding a maximal exact factorization into separable subsystems of potentially multiple qubits apiece. I'll address both questions, because both point to techniques for determining exact subsystem separability.

When we talk about "the states of individual qubits" even when they're potentially entangled, we can mean the "reduced density matrix" representation of the state. For single qubits, it's easy to find this representation: find the expectation values in Pauli Z, X, and Y bases of any individual qubit (even if it's entangled). These are the three components of a representation of the local qubit state on the Bloch sphere. Assuming your original overall quantum state is a "pure state," if the qubit is entangled, then the reduced density matrix Bloch sphere representation will appear as a "mixed state" locally, and the vector will be "off-shell" and interior to the Bloch sphere radius. However, by corollary, if the qubit is separable, the length of the Bloch sphere state vector will be "on-shell," with maximum radius. In classical simulation of quantum mechanics, this is a relatively inexpensive method for checking single-qubit separability (and authors including myself leverage it for a more general approximation technique in arXiv:2304.14969).

As for a maximal factorization of any quantum state, while preserving all entanglement of subsystems, a historically common approach to the problem involves "Schmidt decomposition." The original overall pure state vector is reshaped according to the dimensionality of potential bipartite subsystem boundaries to check, and singular value decomposition (SVD) is performed: if there is only one nonzero singular value (i.e. we say the "Schmidt rank" is 1), then this immediately tells us that the original state can be written as the Kronecker product of two independent subsystems along the tested boundary. According to the usual conventions of singular value decomposition, the first columns of the intermediate matrices called "U" and "V" contain the actual factorized subsystem state vectors. It's computationally "hard," but full factorization of the state could be achieved at least in principle by testing every potential bipartite subsystem boundary, iteratively.

As author of the unitaryfund/qrack quantum computer simulation library and framework, it's worth pointing out that Qrack's "Decompose() method" and "TryDecompose() method" constitute a novel alternative to Schmidt decomposition, to achieve the same result with potentially much less computational space complexity and significantly less time complexity. A notable caveat is that I lack a formal proof of the general correctness of this technique, though it has been used ubiquitously in Qrack internals since 2018 and for over 1.1 million recorded downloads, and no "empirical" evidence has been observed by authors or users of the library that the technique fails in any case that Schmidt decomposition wouldn't. (The algorithm resulted from looking at various examples of states produced via the Kronecker product and attempting to infer a direct inverse in the cases it exists. Formal proof might be trivial, but I honestly simply don't know, and I haven't really worried about it for the years since 2018 that it's never failed in unit tests, integration tests, or any case of field deployment that we needed.)