Expansion of multi-qubit density matrix in the Pauli matrix basis

Question

The single qubit density matrix can be expanded as $$ \rho=\frac{tr(\rho)I+tr(X\rho)X+tr(Y\rho)Y+tr(Z\rho)Z}{2} $$ which can be shown as,

$\rho$ is a positive operator with $tr(\rho)=1$, ie., $\rho=\begin{bmatrix}u&v+iw\\v-iw&1-u\end{bmatrix}$ and $X=\begin{bmatrix}1&0\\0&1\end{bmatrix},Y=\begin{bmatrix}0&-i\\i&0\end{bmatrix},Z=\begin{bmatrix}1&0\\0&-1\end{bmatrix}$ are the Pauli matrices. $$ X\rho=\begin{bmatrix}v-iw&1-u\\u&v+iw\end{bmatrix}\implies tr(X\rho)=2v\\ Y\rho=\begin{bmatrix}-iv-w&-i-iu\\iu&iv-w\end{bmatrix}\implies tr(Y\rho)=-2w\\ Z\rho=\begin{bmatrix}u&v+iw\\-v+iw&-1+u\end{bmatrix}\implies tr(Z\rho)=2u-1\\ 2\rho=\begin{bmatrix}2u&2v+i2w\\2v-i2w&2-2u\end{bmatrix}=\begin{bmatrix}1+(2u-1)&2v+i2w\\2v-i2w&1-(2u-1)\end{bmatrix}\\ =I+2vX+(-2w)Y+(2u-1)Z\\ =tr(\rho)I+tr(X\rho)X+tr(Y\rho)Y+tr(Z\rho)Z\\ \implies \boxed{\rho=\frac{tr(\rho)I+tr(X\rho)X+tr(Y\rho)Y+tr(Z\rho)Z}{2}} $$

How do I prove that the density matrix for $n$ qubits can be written as, $$ \rho=\sum_{\vec{v}}\frac{tr(\sigma_{v_1}\otimes\sigma_{v_2}\otimes\cdots\otimes\sigma_{v_n}\rho)\sigma_{v_1}\otimes\sigma_{v_2}\otimes\cdots\otimes\sigma_{v_n}}{2^n} $$

My Attempt

Lets consider a $2$ qubit case,

We need to prove $\rho=\sum_{i,j} \frac{tr(\sigma_i\otimes\sigma_j\rho)\sigma_i\otimes\sigma_j}{4}$

Lemma: Let $V$ and $W$ be normed spaces. If $\{v_i\}\in V$ and $\{w_j\}\in W$ are linearly independent in $V$ and $W$ respectively, then $\{v_i\otimes w_j\}$ is linearly independent in the algebraic tensor product $V\otimes W.$

Since we can define $\mathbb{C}^{16}$ an isometric to the matrix space $M_{4\times 4}(\mathbb{C})$, lemma can be used to prove $\sigma_i\otimes\sigma_j$ form an orthogonal basis for the $2$ qubit density matrix, i.e., $\rho=\sum_{i,j} c_{i,j}\sigma_i\otimes\sigma_j$

So, how do I prove that the coefficients $c_{i,j}=tr(\sigma_i\otimes\sigma_j\rho)$ ?

Answer 1

"Trace out" $\sigma_i \otimes \sigma_j$

$$\rho = \sum c_{ij}\sigma_i \otimes \sigma_j$$

$$\rho \ \sigma_m \otimes \sigma_n = \sum c_{ij} \ \sigma_i \sigma_m \otimes \sigma_j \sigma_n$$

\begin{align} \mathrm{Tr \{\rho \ \sigma_m \otimes \sigma_n\}} &=\mathrm{Tr}(\sum c_{ij} \ \sigma_i \sigma_m \otimes \sigma_j \sigma_n)\\ &=\sum \mathrm{Tr}(c_{ij} \ \sigma_i \sigma_m \otimes \sigma_j \sigma_n)\\ &=\sum c_{ij} \ \mathrm{Tr}(\sigma_i \sigma_m \otimes \sigma_j \sigma_n)\\ &=\sum c_{ij} \ \mathrm{Tr}(\sigma_i \sigma_m)\mathrm{Tr}(\sigma_j \sigma_n)\\ &= 4\sum c_{ij} \delta_{im} \delta_{jn}\\ &=4c_{mn} \end{align}

$$c_{mn} = \frac{\mathrm{Tr \{\rho \ \sigma_m \otimes \sigma_n\}}}{4}$$

Where I have used the property $\mathrm{Tr}(\sigma_a \sigma_b) = 2\delta_{ab}$, where $\sigma_a$ belongs to the set $\{1_2, \sigma_1, \sigma_2, \sigma_3 \}$, i.e. identity followed by the Pauli matrices.