From any operator you get a matrix when using a basis. If you have a given matrix you need to tell the basis to reconstruct the operator. In QT you have a Hilbert space, and usually you work with a orthonormal basis, ##|u_k \rangle##, ##k \in \mathbb{N}## or also quite often in a "generalized basis" when you describe the Hilbert space in terms of wave functions like in position or momentum representation, but let's stick to the simpler case of a discrete orthonormal basis; for a single particle in non-relativistic QT you can take the harmonic-oscillator energy eigenstates as a nice example.
Now given an operator ##\hat{O}##, you can write it with help of the basis in terms of its "matrix elements",
$$O_{jk}=\langle u_j|\hat{O} u_k \rangle.$$
The orthonormal set is also complete, i.e., you can write the unit operator using the completeness relation:
$$\sum_{k=1}^{\infty} |u_k \rangle \langle u_k \rangle=\hat{1}.$$
Using this completeness relation twice, you can easily reconstruct the operator using the same basis taken to evaluate its matrix elements:
$$\hat{O}=\sum_{j,k=1}^{\infty} |u_j \rangle \langle u_j|\hat{O} u_k \rangle \langle u_k|=\sum_{j,k=1}^{\infty} |u_j \rangle \langle{u}_k| O_{jk}.$$
The most convenient set for such a basis is to choose an eigenbasis of the operator in question. If the operator is self-adjoint (as is the case for the density operator), the corresponding eigenvectors can be made to an orthornomal set. Then you have
$$\hat{O} |u_j \rangle=o_j |u_j \rangle,$$
where ##o_j## is a eigen value of the operator and ##u_j## a correspondin eigenvector. Then the matrix elements simplify very much:
$$o_{jk} = \langle u_j|\hat{O} u_k \rangle = o_k \langle u_j|u_k \rangle=o_k \delta_{jk}.$$
Note that here no Einstein summation convention is applied! This tells you that in the eigenbasis the matrix representing the operator is diagonal with the eigenvalues as entries on the diagonal.
Then one of the sum in the general decomposition formula can be done, and it simplifies to
$$\hat{O}=\sum_{jk=1}^{\infty} |u_j \rangle \langle u_k | o_{jk} = \sum_{jk=1}^{\infty} |u_j \rangle \langle u_k | o_k \delta_{jk}=\sum_{j=1}^{\infty} |u_j \rangle \langle u_j| o_j.$$
All this applies to the density operator (a better name is "statistical operator"). Then the eigen values are by definition non negative and their sum adds to 1, but that doesn't change the above given general framework.