How to efficiently calculate the inverse of a Kronecker product?

Question

This is a follow-up question to a previous question I had, where the correct answer was to use the Kronecker product.

Given, for example, a vector representing two qubits $$\begin{bmatrix}0 \\ 1 \\ 0 \\ 0\end{bmatrix}$$ is there an algorithm to decompose it into two vectors $$\begin{bmatrix}1 \\ 0\end{bmatrix}$$ and $$\begin{bmatrix} 0 \\ 1\end{bmatrix}$$

If so, does this generalize to N qubits? I suspect the procedure doesn't scale well computationally (hence the entire point of quantum computing) but I'm still curious to know whether and what form the algorithm may take.

Answer 1

You specifically ask about qubits, so I'll keep it to that. Imagine you have a state $$ |\psi\rangle=\sum_{x\in\{0,1\}^n}a_x|x\rangle. $$ You can choose to look at each qubit. I'll take the first qubit for the sake of simplicity. We have that $$ |\psi\rangle=|0\rangle\sum_{y\in\{0,1\}^{n-1}}a_{0y}|y\rangle+|1\rangle\sum_{y\in\{0,1\}^{n-1}}a_{1y}|y\rangle $$ which we can equally well think of as $$ |\psi\rangle=c_0|0\rangle|\phi_0\rangle+c_1|1\rangle|\phi_1\rangle, $$ where $\langle\phi_0|\phi_0\rangle=\langle\phi_1|\phi_1\rangle=1$. If it happens that $\langle\phi_0|\phi_1\rangle=e^{i\gamma}$, then this is the same as $$ |\psi\rangle=(c_0|0\rangle+c_1e^{i\gamma}|1\rangle)|\phi_0\rangle, $$ which is now in its separable form, as desired. You keep repeating this qubit-by-qubit. Of course, as the number of qubits, $N$, increases, you're dealing with exponentially many coefficients, so the calculation is not efficient, although it is procedurally straightforward.

In terms of actually performing the computation, it may be easier to rely on the partial trace, which is built in to many packages. You'd calculate, for example, $$ \rho_1=\text{Tr}_{2,3,\ldots N}|\psi\rangle\langle\psi|, $$ the reduced density matrix on the first qubit. If $\text{Tr}\rho_1^2=1$, then this state is separable. Moreover, $\rho_1$ is a rank 1 projector, and the state that it projects onto is the separable decomposition of that qubit.

Answer 2

Consider an arbitrary vector $\boldsymbol v\equiv(v_i)_{i=1}^N\in\mathbb C^{nm}$ of length $n m$ (if you care only about qubits, just fix $n=m=2$).

For $\boldsymbol v$ to have an $(n,m)$ partition means that we can write it as $\boldsymbol v=\boldsymbol u\otimes \boldsymbol w$ for some $\boldsymbol u\in\mathbb C^n$ and $\boldsymbol w\in\mathbb C^m$. If this is the case, then the components of $\boldsymbol v$ must be related to those of $\boldsymbol u,\boldsymbol w$ via $$v_{ij}=u_i w_j. \tag A$$ The first thing to notice is that the tensor product operation $(\boldsymbol u,\boldsymbol w)\mapsto\boldsymbol v$ is not injective. In particular, for any $\lambda\in\mathbb C$ we have $(\lambda \boldsymbol u)\otimes \boldsymbol w=\boldsymbol u\otimes(\lambda \boldsymbol w)$, and thus given a separable $\boldsymbol v$ we can only recover its "reduced vectors" up to their individual norms and phase factors. Because we are interested in quantum states, this means that we don't need to worry about neither the normalisation nor the phase of the vectors $\boldsymbol u,\boldsymbol w$.

A remarkable property of (A) is that, for any pair of indices $j\neq j'$ and $i$, we have $v_{ij}/v_{ij'}=w_j/w_{j'}$, where, you might notice, the $i$ index does not appear on the right-hand side.

This property characterises separable vectors: check that $v_{ij}/v_{ij'}$ does not depend on $i$ for all $j,j'$, and if so, define $w_i\equiv v_{i1}/v_{ij}$, where we are using the fact that $\boldsymbol w$ is defined up to normalisation to assume $w_1=1$ (when some of the elements $v_{ij}$ we might not be able to assume this, but those cases are usually trivial to handle separately). Once you have $\boldsymbol w$, you can find $u_i$ by simply computing $v_{ij}/w_j$.

For example, if $\boldsymbol v=(1,1,2,2)^T$, then $$v_{11}=v_{12}=1, \qquad v_{21}=v_{22}=2,$$ and thus $v_{12}/v_{11}=v_{22}/v_{21}=1$. You can thus define $\boldsymbol w=(1,1)$. You then compute $u_1=v_{11}/w_1=1$ and $u_2=v_{21}/w_1=2$, which tells you that $\boldsymbol u=(1,2)$.

For example, if $\boldsymbol v=(0,1,0,0)^T$ and $n=m=2$, then by definition $$v_{11}=v_{21}=v_{22}=0 \qquad\text{and}\qquad v_{12}=1.$$ Because $v_{11}=0$, we cannot use directly the procedure above, but the fact that $v_{2j}=0$ for all $j$ tells you immediately that $\boldsymbol u=(1,0)$, and then finding $\boldsymbol v$ is easy.

As for how numerically efficient this particular method is, I'm not so sure. The main problem is that all these divisions will be prone to numerical instability for small values of the amplitudes. Methods that are based on computing entropies might work better, but I'm not sure.

Answer 3

Mathematica Code :

Inverting the boolean vector to qubit vector, no checks are added, the input vector should be of the order $2^n$, the resultant vector would be of size $n$

booleanToQubit[v0_]:=Module[{v=v0},
        r={};
        While[
        Length[v]>=2,
        t1=Partition[v,2];
        x={};

        If[Length[Position[t1,{1,0}]]==1,
        t1[[Flatten@Position[t1,{1,0}][[1]]]]=1;
        x={1,0}
        ];  
        If[Length[Position[t1,{0,1}]]==1,
        t1[[Flatten@Position[t1,{0,1}][[1]]]]=1;
        x={0,1}
        ];  

        r=Insert[r,x,1];
        v=t1/.{{0,0}->0};
        ];  
        r/.{{0,1}->0,{1,0}->1}
]

booleanToQubit[{0, 1, 0, 0, 0, 0, 0, 0}]
{1, 1, 0}

booleanToQubit[{0, 1, 0, 0}]
{1, 0} /. {1 -> {1, 0}, 0 -> {0, 1}}, qubit representation
{{1, 0}, {0, 1}}, your input