How to write a two qubit state as "diagonal" in the basis of Pauli matrices?

Any two qubit density matrix can be written as

$$\rho = \frac{1}{4} \sum_{n,m = 0}^{3} M_{nm} (\sigma_n \otimes \sigma_m),$$

where $$\sigma_\mu$$'s are the identity and Pauli matrices.

Is it possible to do a local unitary transformation on the system, like $$U_1 \otimes U_2$$, which will transform $$\rho$$ into $$\rho'$$ which looks like this:

$$\rho' = \frac{1}{4} \left( \sigma_0 \otimes \sigma_0 + M^{'}_{11} \sigma_1 \otimes \sigma_1 + M^{'}_{22} \sigma_2 \otimes \sigma_2 + M^{'}_{33} \sigma_3 \otimes \sigma_3 \right)$$

i.e. the matrix whose elements are $$\{M_{nm}\}$$ is diagonal?

This is indeed possible. Consider the SVD of $$M$$, written as $$M_{nm}=\sum_k s_k \langle n |u_k\rangle \overline{\langle m| v_k\rangle} \equiv \sum_k s_k u_{nk} \bar v_{mk}$$ for some pair of orthonormal bases $$\{|u_k\rangle\}$$ and $$\{|v_k\rangle\}$$ and positive reals $$s_k>0$$. More concisely, this amounts to $$M=\sum_k s_k |u_k\rangle\!\langle v_k|$$. Then $$4\rho = \sum_{knm} s_k u_{nk} \bar v_{mk}(\sigma_n\otimes \sigma_m) = \sum_k s_k \left( \sum_n u_{nk} \sigma_n \right) \otimes \left( \sum_m \bar v_{mk} \sigma_m \right).$$ Now observe that any unitary $$U$$ acts via the adjoint on Pauli matrices (or any other orthonormal basis of operators) as a unitary operator, meaning $$U\sigma_i U^\dagger = \sum_j v_{ji} \sigma_j$$ for some unitary $$\tilde U$$ with components $$v_{ji}$$. With this, we see immediately that we can write $$4\rho = \sum_k s_k (U\sigma_k U^\dagger)\otimes(V\sigma_k V^\dagger)$$ for some pair of unitaries $$U,V$$, and thus we conclude $$4(U\otimes V)^\dagger\rho (U\otimes V)=\sum_k s_k (\sigma_k\otimes\sigma_k).$$
Some (maybe too many?) details about the above relation for $$U\sigma_i U^\dagger$$ can be found here and here. I'm sure I've seen more direct treatments of it on some SE site, but I can't find the post right now.
1. After writing the SVD of $$M$$, I implicitly assumed that the numbers $$u_{nk}$$ are the elements of some unitary matrix $$U\equiv (u_{nk})_{nk}$$. This is not always the case. We can however assume that is without loss of generality, by simply allowing $$s_k=0$$ for some $$k$$. This way we always have $$4$$ elements in the sum, and $$U$$ is always unitary.
2. Any special unitary $$U\in\mathbf{SU}(2)$$ can be written as a sum of Pauli operators as $$U = \alpha_0 I + i \sum_{j=1}^3 \alpha_j \sigma_j,$$ with $$\alpha_i\in\mathbb{R}$$ and $$\sum_{i=0}^4\alpha_i^2=1$$. From this decomposition, you can compute explicitly $$U \sigma_i U^\dagger = \left(\alpha_0^2-\sum_{j=1}^3\alpha_j^2\right)\sigma_i + 2 \alpha_i (\boldsymbol\alpha\cdot\boldsymbol\sigma) +2i\alpha_0 (\boldsymbol\sigma\times\boldsymbol\alpha)_i.$$