Extracting $\mathbf b$ from $M = a I + \mathbf b \cdot \mathbf S$, when $S_i$ are higher spin matrices? This is a cross-post of a question that I posted on the Math SE, that did not get any answers there. It is fundamentally a mathematics question, but it pertains to spin matrices, which many Physics SE users have a lot of knowledge about, so I thought I might try here. My end use case is also in physics, rather than pure mathematics.
The original question
Any $2 \times 2$ hermitian matrix $M$ can be written
$$M = a I + \mathbf b \cdot \boldsymbol \sigma,\tag{1}$$
where $a \in \mathbb R$, $\mathbf b \in \mathbb R^3$, $\boldsymbol \sigma$ is the Pauli vector, and $\mathbf b \cdot \boldsymbol \sigma := \sum b_i \sigma_i$. This is because the identity matrix and the Pauli matrices make up a basis for the real vector space of $2 \times 2$ hermitian matrices. Moreover, given $M$, one can extract $a$ and $\mathbf b$ via
$$a = \frac{1}{2} \mathrm{tr}(M), \quad b_i = \frac{1}{2} \mathrm{tr}(\sigma_i M)\tag{2}$$
because the Pauli matrices are traceless and obey the identity
$$\sigma_j \sigma_k = \delta_{jk} I + i \epsilon_{jkl} \sigma_l,\tag{3}$$
where we use implicit summation over the repeated index $l$.
A generalization of the Pauli matrices that is of particular interest to physicists is the higher spin matrices. Hence, say we replace the Pauli vector in (1) with a "vector" $\mathbf S$ of $d \times d$ spin matrices, yielding
$$M = a I + \mathbf b \cdot \mathbf S,\tag{4}$$
where now, of course, $M$ and $I$ are $d \times d$. (This is clearly not a general $d \times d$ hermitian matrix for $d > 2$, but more restricted.) Is there then some simple way, similar to (2), to extract $\mathbf b$ given $M$?
I tried the cases $d = 3, 4$ in Mathematica and found that, indeed,
$$b_i \propto \; \mathrm{tr}(S_i M).\tag{5}$$
So I suspect this works in general, but my initial attempts at a simple proof have failed. The identity (3) does not hold for the higher spin matrices, because, while they do obey
$$[S_j, S_k] = i \epsilon_{jkl} S_l\tag{6}$$
much like the Pauli matrices (up to normalization), they do not obey the same anticommutation relations.
 A: Because the higher spin matrices generate a representation of $SO(3)$, we have :
$$\operatorname{Tr}(S_iS_j) \propto \delta_{ij} \tag{1}$$
(because $\delta_{ij}$ is the only invariant rank 2 tensor)
To find the normalization constant in $(1)$, contract $i$ and $j$ and express the RHS in terms of the total spin.
On top of that, we know that $\operatorname{Tr}(S_i) = 0$. Therefore :
$$\operatorname{Tr}(S_i(aI + B_jS_j)) = B_j\operatorname{Tr}(S_iS_j)\propto B_i$$
Edit
Equation $(1)$ is a classical result from the theory of simple Lie algebras (that the Killing form is positive definite (or negative definite if you choose anti-hermitian generators as mathematicians do)). Therefore you can choose generators so that $(1)$ is satisfied.
Now, in a given representation $R$, let us write $\newcommand{\Tr}{\operatorname{Tr}} \Tr_R(S_iS_j) = T(R)\delta_{ij}$ where $T$ is called the intex of the representation. $\mathbf S^2 = S_iS_i = C(R)$ is the quadratic Casimir (which I called total spin above).
Contracting $i$ and $j$ above we get :
$$\dim(R)C(R) = T(R)\dim(G)$$
In the case of $SU(2)$ or $SO(3)$, the spin $J$ representation has $\dim(R) = 2J+1$, $C(G) = J(J+1)$ so we get :
$$T(R) = \frac{J(J+1)(2J+1)}{3}$$
