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)$ ?