Why should we diagonalize the density matrix? Why is the diagonal representation more special? Suppose I have
$$\rho=\sum p_k \vert \phi_k \rangle \langle \phi_k\vert = \sum \lambda_k \vert \lambda_k \rangle \langle \lambda_k\vert$$
where I call the second decomposition "orthonormal", meaning $\langle \lambda_l\vert \lambda_m\rangle =\delta_{lm}$, and the numbers before the projectors are the probabilities.
Why is the orthonormal decomposition more special than a generic one? By "special" I am referring to, for example, the fact that the Von Neumann entropy $S$ of the density operator and the Shannon entropy $H$ of the state distribution coincide only in the orthonormal decomposition  $$H(\lbrace p_i \rbrace)=\sum p_i \log \frac{1}{p_i} \geq H(\lbrace \lambda_i \rbrace) = S(\rho)=-\operatorname{tr}\rho \log \rho.$$
Can it be said that othornomal states are more "independent" than those which are not? Can this be made mathematically precise?
 A: Well, this seems to be a linear algebra question. 
Firstly, mathematically, when you want to calculate the trace, i.e.
\begin{align}
tr(\cdot) = \sum_{\lambda}\langle\lambda|\cdot|\lambda\rangle
\end{align}
This calculation of course requires the $|\lambda\rangle$-basis is an orthogonal one -- think about in the original matrix language in math, any matrix, when assigned with a linear space, implicitly uses a orthogonal basis. 
As you have already mentioned, the definition of entanglement entropy uses trace:
\begin{align}
S &= -tr\big[\rho\log{\rho}\big] \\
&= -\sum_{\lambda}\langle\lambda|\rho\log{\rho}|\lambda\rangle
\end{align}
Firstly, for sure you could use any basis $|\phi\rangle$, orthogonal or un-orthogonal, to calculate $S$, as long as you can simply calculate each:
\begin{align}
\langle\lambda|\phi\rangle
\end{align}
And there are more than one orthogonal-basis, for sure, related by orthogonal transformations. But if you are not using the "eigen-basis" in which $\rho$ has a diagonal form, then you'll have to calculate the above quantity anyway -- that means difficulty of calculation has not been reduced at all. But, say, if you could use the diagonal basis, then the trace calculation would be simply a summation of diagonal elements of $\rho\log{\rho}$.
A: You have already discussed the convenience of the diagonal representation for entropy calculations. I will add to that from a state preparation point of view to motivate the diagonal representation. Say you're in a lab and you want to make a mixed state $\rho$. The ensemble ambiguity theorem of quantum mechanics tells you that 
1) There are infinitely many ways of preparing this state.
2) No matter how you prepare the state, no experiment can distinguish between the methods of preparation.
One of the simplest choices for state preparation you could make is roll a dice (with the appropriate probability distribution) and choose a pure state from an orthogonal set depending on the dice outcome. But what you have done here is write $\rho$ as a convex combination of only orthogonal pure states i.e. $\rho = \sum_i p_{i}\vert i\rangle\langle i\vert$. This method of preparation requires the diagonal representation. 
Moreover, you prepare a pure state from some orthonormal set each time with only classical randomness from the dice that decides which state. Hence, all the entropy arguments that connect the von Neumann and Shannon entropies that were detailed in the other answer immediately follow. 
