How can projection operators be expressed in form $\frac{1}{d} (I + \sum_i r_i \lambda_i)$? How can projection operator be expressed in form $\frac{1}{d} (I + \sum_i r_i \lambda_i)$?
I was reading a paper and found out that the density matrix in $d$-dimensional Hilbert Space can be expressed in form
$$\frac{1}{d} (I + \sum_i r_i \lambda_i),$$
where $i = 1,2,\ldots,d^2-1$, $r_i$ are reals (with certain conditions), $\lambda_i$ is the generators of $SU(d)$ algebra, they are orthogonal in the sense that
$$tr(\lambda_i\lambda_j) = \delta_{ij}.$$
I was wondering is there a way to explicitly choose $n$ sets of $\{r_i\}$ to form a complete set of orthogonal basis, i.e. projection operators, just like in $2$-dimensional case, where
$$\frac{1}{2} (I \pm (a_1\sigma_1 + a_2 \sigma_2 + a_3 \sigma_3))$$
are the two orthogonal  projection operators, with $(a_1,a_2,a_3)$ on unit sphere.
 A: The generators you mention form a complete basis for traceless hermitian matrices. Thus, to express a hermitian matrix that has a non zero trace in this basis, you have to augment the basis with one more matrix that is trace orthogonal to the others, and the identity does the trick. You can verify that the number of matrices in your basis agrees with the number of degrees of freedom of a hermitian matrix.
Another useful fact is that the generators of $U(N)$ are exactly the matrices you listed. The identity matrix take care of the phase of unitary matrices.
Btw, the "dot product" for this basis is just the trace of the product of the two matrices, so you can use this to decompose an arbitrary Hermitian matrix in your basis.
A: Since I started writing it, I will complete this answer to add some extra mathematical detail. As stated by lionelbrits, the $d^2-1$ Hermitian generators of $\mathrm{SU}(d)$ are traceless and orthonormal with respect to the Hilbert-Schmidt inner product
$$(A,B) = \mathrm{Tr} (A^{\dagger}B).$$
Therefore, along with the identity, these operators form a complete, orthonormal basis in the $d^2$-dimensional Hilbert space of operators acting on a $d$-dimensional Hilbert space of vectors. Any operator can be expanded as an appropriate linear combination of the basis operators; a Hermitian operator will have real expansion coefficients. 
It is convenient to normalise the basis so that
$\mathrm{Tr} (\lambda_i\lambda_j) = d\delta_{ij}.$
This is achieved by multiplying by a factor $\sqrt{d}$ the generators normally used by mathematicians. This convention is consistent with the OP's use of Pauli matrices in the case $d=2$. With this normalisation, any rank-one projector $\Pi = \lvert \psi\rangle\langle\psi\rvert$ can be expanded as
$$\Pi  = \frac{1}{d}\left( I + \sum_{i=1}^{d^2-1}r_i\lambda_i \right),$$
where $r_i = \mathrm{Tr}(\Pi \lambda_i)=\langle\psi\rvert\lambda_i\lvert\psi\rangle$. The condition $\mathrm{Tr} (\Pi^2) = 1$ leads to the restriction
$$ \sum_{i=1}^{d^2-1}r_i^2 = d-1,$$
meaning that the $r_i$ are coordinates on a $(d^2-2)$ -sphere of radius $\sqrt{d-1}$.
If two rank-one projectors with expansion coefficients $\{r_i\}$ and  $\{s_i\}$ are orthogonal, one finds the condition $\sum_i r_is_i = \mathbf{r}\cdot\mathbf{s} = -1$. When $d=2$, this implies that $\mathbf{r} = -\mathbf{s}$, i.e. the projectors are described by antipodal points on the unit sphere. For $d>2$ this is no longer the case: instead one finds that the angle between $\mathbf{r}$ and $\mathbf{s}$ is
$$\theta = \arccos \left(\frac{-1}{d-1}\right)$$
if the projectors are orthogonal. This illustrates that the nice geometrical intuition afforded by the Bloch sphere when $d=2$ does not carry over so simply to higher dimensions.
