Consider an elliptic self-adjoint pseudodifferential operator $A$ acting on $m$-columns of half-densities on a closed manifold $M$, whose principal symbol is assumed to have simple eigenvalues. Relying on a basis of pseudodifferential projections commuting with $A$, we construct an almost-unitary pseudodifferential operator that diagonalizes $A$ modulo an infinitely smoothing operator. We provide an invariant algorithm for the computation of its full symbol, as well as an explicit closed formula for its subprincipal symbol. Finally, we give a quantitative description of the relation between the spectrum of $A$ and the spectrum of its approximate diagonalization, and discuss the implications at the level of spectral asymptotics.