Prove that the ensemble decompositions of a density operator $\rho$ are related by a unitary operator

Question

In my class it was told that ensemble decompositions of a density operator $\rho$ are not unique, but that the ones that exist are related by a unitary operator. I'm trying to prove this, but I get stuck somewhere along the way.

Lets begin by assuming two different decompositions of density operator $\rho$: $\rho = \sum_{j=1}^n{p_j|\psi_j\rangle\langle\psi_j|} = \sum_{k=1}^m{q_k|\phi_k\rangle\langle\phi_k|}$

Now, these two decompositions live in a Hilbert space $\mathcal A$. We can then define a purification of both, using a system described by a Hilbert space $\mathcal B$ of dimension $k=\max (n,m)$, so that we get $|\Psi_1\rangle_{\mathcal A\mathcal B} = \sum_{j=1}^n{\sqrt{p_j}|\psi_j\rangle \otimes |b_j\rangle}$ and $|\Psi_2\rangle_{\mathcal A\mathcal B} = \sum_{k=1}^m{\sqrt{q_k}|\phi_k\rangle \otimes |b_k\rangle}$.

Now, here we can use that as these pure states are purifications of the same density operator, there must be a unitary $U$ connecting them: $(1_A \otimes U_B)|\Psi_1\rangle_{\mathcal A\mathcal B} = |\Psi_2\rangle_{\mathcal A\mathcal B}$.

Here is where I get stuck. I should be able to use this to prove the unitary relation between the $\psi$ and the $\phi$, but it is not obvious to me how I should do this.

Update: after reviewing the comments to the first question, I should have written that the $\psi$ and $\phi$ states do NOT have to be orthonormal, per se.

Answer 1

Result for one orthogonal and one non-orthogonal decomposition

We will prove the following Theorem:

Let $\rho=\sum_{i=1}^Np_i\lvert\phi_i\rangle\langle\phi_i\rvert$ be a eigenvalue decomposition and $\sigma=\sum_{i=1}^Mq_i\lvert\psi_i\rangle\langle\psi_i\rvert$ ($M\ge N$). Then, $\rho=\sigma$ if and only if \begin{equation} \sqrt{q_j}\lvert\psi_j\rangle=\sum_i v_{ij}\sqrt{p_i}\lvert\phi_i\rangle\ , \end{equation} with $\sum_j v_{ij}v_{i'j}^*=\sum_k \delta_{ii'}$, i.e., $V\equiv(v_{ij})$ is an isometry.

Proof:

The "if" direction is straightforward: \begin{align} \rho&=\sum_{j}q_j\lvert\psi_j\rangle\langle\psi_j\rvert\\ &=\sum_{i,i',j} v_{ij} v_{i'j}^* \sqrt{p_i}\sqrt{p_{i'}} \lvert\phi_i\rangle \langle\phi_{i'}\rvert \\ &=\sum_i p_i\lvert\phi_i\rangle \langle\phi_{i}\rvert = \sigma\ , \end{align} where in the last step we have used that $\sum_j v_{ij}v_{i'j}^*=\delta_{ii'}$.

To prove the converse, let $$ v_{ij} := \langle\phi_i\lvert\psi_j\rangle\,\sqrt{q_j/p_i}\ .$$ Then, $$ \sum_i v_{ij} \sqrt{p_i}\lvert\phi_i\rangle = \sum_i \lvert\phi_i\rangle\langle \phi_i\lvert\psi_j\rangle\sqrt{q_j} = \sqrt{q_j}\lvert\psi_j\rangle\ ,$$ i.e., $v_{ij}$ is the desired basis transformation. Further, $$ \sum_{ii'}\underbrace{\delta_{ii'}\sqrt{p_ip_{i'}}}_{=:a_{ii'}}\lvert\phi_i\rangle\langle\phi_i\rvert = \sum_ip_i\lvert\phi_i\rangle\langle\phi_i\rvert = \sum_jq_j\lvert\psi_j\rangle\langle\psi_j\rvert = \sum_{ii'}\underbrace{\sum_jv_{ij}v^*_{i'j}\sqrt{p_ip_i'}}_{=:b_{ii'}}\lvert\phi_i\rangle\langle\phi_{i'}\rvert\ . $$ Now, since the $\lvert\phi_i\rangle$ are orthogonal (as they form an eigenbasis), the $\lvert\phi_i\rangle\langle\phi_{i'}\rvert$ are linearly independent, and thus, $a_{ii'}=b_{ii'}$, which implies $\sum_jv_{ij}v^*_{i'j}=\delta_{ii'}$.

---

Extension to two non-orthogonal decompositions

If the $\lvert\phi_i\rangle$ do not form an orthonormal basis, we can generalize the theorem by going through an orthonormal basis $\rho=\sum r_k\lvert\chi_k\rangle\langle\chi_k\rvert$: Then, \begin{align} \sqrt{p_i}\lvert\phi_i\rangle&=\sum_k u_{ki}\sqrt{r_k}\lvert\chi_k\rangle\\ \sqrt{q_j}\lvert\psi_j\rangle&=\sum_k w_{kj}\sqrt{r_k}\lvert\chi_k\rangle \end{align} with $\sum_i u_{ki}u_{k'i}^*=\delta_{kk'}$ and $\sum_j w_{kj}w_{k'j}=\delta_{kk'}$. The second equation then yields $\sqrt{r_k}\lvert\chi_k\rangle= \sum_j w_{kj}^*\sqrt{q_j}\lvert\psi_j\rangle$. After inserting this in the first equation, we obtain $$ \sqrt{q_j}\lvert\psi_j\rangle=\sum_i v_{ij}\sqrt{p_i}\lvert\phi_i\rangle $$ with $v_{ij} = \sum_k u_{ki}w_{kj}^*$, i.e., $V=UW^\dagger$ is a partial isometry.

Answer 2

A slightly different spin on the characterisation of the set of ensemble decompositions of a given state.

Let $\rho$ be a generic positive semidefinite operator. I'm not going to assume normalisation here, as it's not necessary for the argument. Being it positive semidefinite, $\rho$ admits an eigendecomposition $$\rho=\sum_k p_k \mathbb{P}(\psi_k),$$ where $p_k>0$ and $\mathbb{P}(\psi_k)\equiv|\psi_k\rangle\!\langle\psi_k|$ are the rank-1 projections on the states $|\psi_k\rangle$. Note in particular that by definition of eigendecomposition we have $\langle\psi_j,\psi_k\rangle=\delta_{jk}$. In general we can have $\sum_k p_k\neq1$ if $\rho$ is not normalised. If $\operatorname{tr}(\rho)=1$ then $\sum_k p_k=1$ and this equation expresses $\rho$ as a convex combination of rank-1 operators. More generally, this equation expresses $\rho$ as an affine combination of such operators. In the following, just replace all mentions of "affine" with "convex" to specialise to the $\operatorname{tr}(\rho)=1$ case.

The goal is to characterise all possible ways to decompose $\rho$ as $\rho=\sum_k \lambda_k \mathbb{P}(\phi_k)$ for some set of $\lambda_k>0$ and not-necessarily-orthogonal states $|\phi_k\rangle$. The trick is to observe that the non-negativity of the coefficients allows to write $\rho=\sum_k \mathbb{P}(\sqrt{\lambda_k}\phi_k)$, which can in turn be written concisely as $\rho= AA^\dagger$ with $$A \equiv \sum_k \sqrt{\lambda_k} |\phi_k\rangle\!\langle k|.$$ In words, $A$ is a linear operator whose columns are the vectors $\sqrt{\lambda_k}|\phi_k\rangle$.

Let me stress that not only any decomposition of $\rho$ in terms of rank-1 operators has this form, but also any linear operator $A$ such that $\rho=AA^\dagger$ corresponds to one such decomposition. In other words, the problem we're discussing is mathematically equivalent to that of characterising the set of linear operators $A$ such that $\rho=AA^\dagger$, for a given $\rho\ge0$. Note that we also don't need $A$ to have the same dimension of $\rho$. If $\rho$ is $d$-dimensional, meaning $\rho\in\operatorname{Lin}(\mathbb{C}^d)$, we can have $A$ of the form $A\in\operatorname{Lin}(\mathbb{C}^n,\mathbb{C}^d)$ for any $n\ge \operatorname{rank}(\rho)$. The idea is here that the dimension of the domain of $A$ corresponds to the number of terms in the decomposition, which can't be smaller than the rank of $\rho$, but can be arbitrarily large (in fact, it can infinite, though I'm only considering finite-dimensional cases here).

So how do we find the linear operators $A$ such that $\rho=AA^\dagger$? The answer is easy thinking in terms of SVDs. This relation imposes $A$ to have singular values equal $\sqrt{p_k}$, and left principal components equal to the eigenvectors of $\rho$. Or more explicitly, the operators we are looking for are all and only those of the form $$A = \sum_k \sqrt{p_k} |\psi_k\rangle\!\langle u_k|,$$ where $|\psi_k\rangle$ are again the eigenvectors of $\rho$. This sum thus clearly contains a number of elements equal to the rank of $\rho$. Furthermore, $\{|u_k\rangle\}$ is here an orthonormal basis for the domain space of $A$, and any such orthonormal basis gives a valid $A$. And any decomposition of $\rho$ corresponds to such a choice of basis. This freedom can be expressed concisely writing that $$A= \sqrt\rho V^\dagger,$$ where we defined the isometry $V\equiv \sum_k |u_k\rangle\!\langle \psi_k|$. Note that this is an isometry, and not in general a unitary operator, because $|u_k\rangle$ can live in a space of higher dimension than $|\psi_k\rangle$.

In conclusion, the set of decompositions for $\rho$ corresponds to the set of linear operators of the form $A=\sqrt\rho V^\dagger$ where $V$ spans all possible isometries. The elements in the decomposition of $\rho$ correspond to the columns of $A$. More explicitly this means that $\rho$ decomposes as $\rho=\sum_k a_k a_k^\dagger$ where $$a_k \equiv A|k\rangle = \sum_j \sqrt{p_j}|\psi_j\rangle \overline{\langle k|V|\psi_j\rangle} = \sum_j \sqrt{p_j} |\psi_j\rangle \bar u_{kj},$$ defining for convenience $U$ as the isometry such that $\langle k|U|j\rangle=\langle k|V|\psi_j\rangle$. This is the formula one usually finds written explicitly.

Answer 3

You certainly want to have $\langle\psi_i|\psi_j\rangle=\delta_{ij}$ and $\langle\phi_k|\phi_l\rangle=\delta_{kl}$, or some equivalent condition ensuring that both ensemble decompositions are optimal. Otherwise, it is easy to come up with many different totally unrelated decompositions.

Because of the optimality condition, both decompositions are eigendecompositions of $\rho$. Hence $n=m$ and $p_i=q_i$, if we order $p_j$ and $q_k$ according to decreasing magnitude. The $|\psi_j\rangle$ and $|\phi_k\rangle$ for a specific eigenvalue are an orthonormal basis basis for the corresponding eigenspace, hence they are related by an unitary operator. For the global unitary operator, extend $|\psi_j\rangle$ and $|\phi_k\rangle$ to an orthonormal basis, and then just map $|\psi_i\rangle$ to $|\phi_i\rangle$. This works even if the $p_i$ and $q_k$ are not ordered, just because of $\langle\psi_i|\psi_j\rangle=\delta_{ij}$ and $\langle\phi_k|\phi_l\rangle=\delta_{kl}$.

So maybe the real task here is to show that being an optimal decomposition implies $\langle\psi_i|\psi_j\rangle=\delta_{ij}$. This is related to properties of the singular value decomposition, which gives a succinct description of the optimal approximations with respect to the Frobenius norm and the spectral norm. Note: The condition $\langle\psi_j|\psi_j\rangle=1$ and $\sum_{j=1}^n{p_j}=1$ is not sufficient to describe this sort of optimality, even so I initially claimed this. So this question turn out to be a bit advanced, because describing the sense in which the decomposition is optimal is non-trivial, if you don't just reference the singular value decomposition for that part.

---

Here is a simple counterexample to the comment by Norbert Schuch that "two decompositions are related by isometries is true regardless of whether the vectors are orthonormal": $\begin{bmatrix}1 & 0 \\ 0 & \epsilon^2\end{bmatrix} = \frac{1}{2}\begin{bmatrix}1 \\ \epsilon\end{bmatrix}\begin{bmatrix}1 & \epsilon\end{bmatrix} + \frac{1}{2}\begin{bmatrix}1 \\ -\epsilon\end{bmatrix}\begin{bmatrix}1 & -\epsilon\end{bmatrix}=\begin{bmatrix}1 \\ 0\end{bmatrix}\begin{bmatrix}1 & 0\end{bmatrix}+\epsilon^2\begin{bmatrix}0 \\ 1\end{bmatrix}\begin{bmatrix}0 & 1\end{bmatrix}$

No isometry can map $(\begin{bmatrix}1 \\ \epsilon\end{bmatrix},\begin{bmatrix}1 \\ -\epsilon\end{bmatrix})$ to $(\begin{bmatrix}1 \\ 0\end{bmatrix},\begin{bmatrix}0 \\ 1\end{bmatrix})$, not even with rescaling, because isometries preserve angles. The first pair of vectors is nearly parallel, while the second pair of vectors is orthogonal.