Representation theorem for matrices (reference request) Motivation. If $A \in \mathbb{C}^{n \times n}$ is self-adjoint (or, more generally, normal), then we all know that
$$
A = \sum_{k=1}^n \lambda_k \, h_k \otimes h_k,
$$
where $\lambda_1,\dots,\lambda_n$ are the eigenvalues of $A$ (counting multiplicities), $(h_1,\dots,h_n)$ is a corresponding orthonormal basis of eigenvectors, and $h_k \otimes h_k$ denotes the matrix given by $(h_k \otimes h_k)x = \langle h_k, x\rangle h_k$ for each $x \in \mathbb{C}^n$ (here I used the "physical" convention that the inner product is linear in the second component).
A representation result. The following result for general (i.e. also non-normal) matrices is - very loosely - reminiscent of the above quoted spectral theorem:
Let $S$ denote the (Euclidean) unit sphere in $\mathbb{C}^n$ and let $\lambda$ denote the surface measure on $S$ (more precisely, we identify $\mathbb{C}^n$ with $\mathbb{R}^{2n}$, consider the surface measure on the unit sphere there and pull it back to $S$). Now, normalize $\lambda$ such that $\lambda(S) = n$.
Theorem. For every matrix $A \in \mathbb{C}^{n\times n}$ we have
$$
A = (n+1) \int_{S} \langle h, Ah \rangle \; h \otimes h \; d \lambda(h) - \operatorname{tr}(A) \, I;
$$
here, $\operatorname{tr}(A)$ denotes the trace of $A$ and $I \in \mathbb{C}^{n\times n}$ denotes the identity matrix.
One can prove the above theorem by using the answers to this MathOverflow question.
The question (a reference request). I have no idea whether the above representation theorem is of any use - but given its very symmetric and rather simple nature, it is natural to suspect that the theorem should already be somewhere in the literature.
So my question is: Do you know any reference where the above representation theorem is stated and proved?
 A: My comments converted to an answer:

1st comment. I know no reference but the proof need not be [as per OP] lengthy — one just computes the two scalars by which the map $\mathscr I:$
$$
\textstyle A\mapsto\mathscr I(A)=\int_S h\langle h,Ah\rangle\langle h,\cdot\rangle d\lambda(h)
\tag1
$$
  acts on the irreducible components of $\mathfrak{gl}(n,\mathbf C)=\mathfrak{sl}(n,\mathbf C)\oplus\mathbf C I$ (Schur’s lemma).

(In detail: write $r$ and $s$ for the scalars in question. Taking the trace of $sI=\mathscr I(I)$ in (1) gives $s=1$. It follows that we have
$
\mathscr I(A)= r\left(A-\tfrac{\operatorname{tr}A}{n}I\right)+\tfrac{\operatorname{tr}A}{n}I
$
and hence
$$
\operatorname{tr}(\mathscr I(A)B)=r\operatorname{tr}(AB)+\tfrac{1-r}n\operatorname{tr}(A)\operatorname{tr}(B)\rlap{\qquad\quad\forall A, B.}
\tag2
$$
Writing this out for $(A,B)=(E_{12},E_{21})$, resp. $(E_{11},E_{22})$ where $E_{ij}=e_i\langle e_j,\cdot\rangle$, one obtains that $\int_S|\langle e_1,h\rangle|^2|\langle e_2,h\rangle|^2d\lambda(h)$ equals both $r$ and $\frac{1-r}n$. Therefore $r=\frac1{n+1}$, and with that (2) becomes (3) below. QED)

2nd comment. For a reference: your desired formula is equivalent to
  $$
\mathrm{tr}(AB)+\mathrm{tr}(A)\mathrm{tr}(B)=(n+1)\int_S\langle h,Ah\rangle\langle h,Bh\rangle\,d\lambda(h)
\tag3
$$
  which is e.g. (3.8) of Gibbons (1992), taking differences of notation into account. (The integrand descends to $P^{n-1}(\mathbf C)$ where he uses the measure of total volume $\pi^{n-1}/(n-1)!$) I think it should also follow from Archimedes-Duistermaat-Heckman (1982, Prop. 3.2).

More context: Generally I think you’ll find many similar formulas in the (somewhat repetitive) literature on “coherent states” or “quantum mechanics as classical mechanics on $P\mathscr H$”. Also compare Schur’s proof of his orthogonality relations (1924, p.199 or Bröcker-tom Dieck 4.5i), specialized to the adjoint representation.
A: This is an easy consequence of the $k=2$ case of the complex version of the Isserlis-Wick theorem for moments of Gaussian measures, i.e., the identity
$$
\int_{\mathbb{C}^n} z_{i_1}\cdots z_{i_k}\ {\bar{z}}_{j_1}\cdots {\bar{z}}_{j_k}
\ e^{-|z|^2}
\prod_{a=1}^{n}\frac{d(\Re z_a) d(\Im z_a)}{\pi}\ 
=\ \sum_{\sigma\in\mathfrak{S}_k} \delta_{i_1 j_{\sigma(1)}}\cdots \delta_{i_k j_{\sigma(k)}}\ .
$$
Going to spherical coordinates produces the integral $\int_{S}\cdots d\lambda(h)$, while the sum over the permutation $\sigma$ gives the other two terms of the wanted identity.
Indeed, it is good to see this really as an integral over $\mathbb{C}\mathbb{P}^{n-1}$
for the Fubini-Study metric, but appealing to Duistermaat-Heckman is not necessary (as known to François).
A: Here is a down-to-earth proof. Remark that for complex matrices, proving $A=B$ is equivalent to proving $\langle x,Ax\rangle=\langle x,Bx\rangle$ ; the superiority of the complex numbers over the real ones !
We must prove that
$$\langle x,Ax\rangle=(n+1)\int_S|\langle h,x\rangle|^2\langle h,Ah\rangle d\lambda(h)-({\rm Tr}\,A)|x|^2.$$
By rotational invariance, it is enough to prove this for $x=\vec e_1$, that is
$$a_{11}=(n+1)\int_S|h_1|^2\langle h,Ah\rangle d\lambda(h)-{\rm Tr}\,A.$$
This amounts to verifying the following identities, all of which being classical:
$$\int_S|h_1|^4d\lambda(h)=\frac2{n+1},\quad\int_S|h_1|^2|h_2|^2d\lambda(h)=\frac1{n+1}$$
and
$$\int_S|h_1|^2h_j\bar h_kd\lambda(h)=0,\quad j\ne k.$$
