We will use $D^l$ for irreducible representations and $V_l$ for the corresponding vector space.
Given two irreducible representations of $SO(3)$, $D^{l_1}$ and $D^{l_2}$ ($l_1$ and $l_2$ are the same as $j_1$ and $j_2$ in Clebsch-Gordan coefficients). Their tensor product $D^{l_1} \otimes D^{l_2}$ may be a reducible representation of $SO(3)$. Indeed, we have
$$D^{l_1} \otimes D^{l_2} = \bigoplus_{| l_1 – l_2 | \leq l \leq l_1 + l_2} D^l$$
The identity is basically representing the same vector space in two different ways. The left hand side corresponds with the natural choice of basis of $V_{l_1} \otimes V_{l_2}$ induced by that of $V_{l_1}$ and $V_{l_2}$. The right hand side is a direct sum of subspaces $V_l$ of $V_{l_1} \otimes V_{l_2}$, each is invariant under $(D^{l_1} \otimes D^{l_2})(g)$ for any $g \in SO(3)$.
Clebsch-Gordan coefficients are basically a change of basis in the sense that they form a transition matrix that changes from the natural basis of $V_{l_1} \otimes V_{l_2}$ to the basis of $V_l$'s.
Now we still need to choose a basis for each $V_l$. My question is: what basis is chosen for each $V_l$ in the context of Clebsch-Gordan coefficients and why are they chosen this way or are they considered natural?
Best Answer
To fix basic Lie group notation: $D^l= \exp (i\theta ~\hat {\mathbf n}\cdot {\mathbf L})$ is a rotation matrix, where the three matrix generators L obey the Lie algebra so(3); and, ipso facto, commute with the quadratic Casimir invariant $\mathbf {L\cdot L}$, whose eigenvalues are all $l(l+1)$ for integer or half-integer $l$, dubbed "spin", for an irreducible representation. (This invariant is proportional to the identity matrix for an irrep.)
Your basic matrix reduction $$D^{l_1} \otimes D^{l_2} = \bigoplus_{| l_1 - l_2 | \leq l \leq l_1 + l_2} D^l$$ consists of $(2l_1+1)(2l_2+1)$-by-$(2l_1+1)(2l_2+1)$ matrices on the left-hand-side acting on the uncoupled basis for the vector space $V_{l_1} \otimes V_{l_2}$ of dimensionality $(2l_1+1)(2l_2+1)$.
The right-hand-side acts on the reduced vector space $\bigoplus_{| l_1 - l_2 | \leq l \leq l_1 + l_2} V_l$ of the same dimensionality, where each $D^l$ acts nontrivially only on the subspace $V_l$, of dimensionality $2l+1$, with null action on all other subspaces. The quadratic Casimir invariant ${\mathbf L}\cdot {\mathbf L}$ has eigenvalue $l(l+1)$ on each different subspace $V_l$, where $l$ varies between the limits indicated in the range of the direct sum.
So, this spherical basis is chosen by the quadratic Casimir invariant, and is considered natural, since it represents disjoint action of the above "huge" representation matrix you wrote: the right-hand-side huge matrix is a sparse collection of smaller block matrices on the diagonal, each of size $2l+1$-by-$2l+1$, for each of the $l$'s specified by the limits given. It represents overall rotation properties of the coupled composite system involved, in which you are presumably interested. You are interested in the properties of the box, not its contents.
The Clebsch-Gordan orthogonal matrix converts the l.h.s. to the r.h.s. basis outlined, where action of the Ls is straightforward for routine computations.
You may illustrate this for $l_1=l_2=1/2$, composition of two doublets, to yield a singlet, $l=0$, directly summed with a triplet, $l=1$.
SU(2) double covers SO(3), but they share the same lie algebras, su(2)=so(3), so all your physics applications would be identical; further discussion of topology would disorient you rather than help you; recall your states are labelled by their Lie algebra eigenvalues l. It is vastly more important to actually learn how to deal with the concrete matrices as exemplified in the suggested illustration above! These are everyday-life issues.
Response to comment hint request:
Typically, you want to find the coproduct generators of your spin-1 $\otimes$ spin-1 generators, $$J_1=\frac{1}{\sqrt 2}\left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right)$$ $$J_2=\frac{i}{\sqrt 2}\left( \begin{array}{ccc} 0 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & 0 \end{array} \right)$$ $$J_3= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{array} \right),$$ with Casimir ${\mathbf J}\cdot {\mathbf J}=2 1\!\!1$.
That is, calculate the 9×9 matrices $$ \Delta (J_i)= J_i\otimes {\mathbb I} + {\mathbb I} \otimes J_i\otimes , $$ which evidently satisfy the same su(2) Lie algebra--note how they relate to the direct product of group matches you inefficiency focussed on, first.
So, how do you reduce these to $$ J_i \oplus S_i \oplus 0, $$ where the Ss are the spin 2 generators $$S_1=\left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & \sqrt{\frac{3}{2}} & 0 & 0 \\ 0 & \sqrt{\frac{3}{2}} & 0 & \sqrt{\frac{3}{2}} & 0 \\ 0 & 0 & \sqrt{\frac{3}{2}} & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ \end{array} \right)$$ $$S_2=\left( \begin{array}{ccccc} 0 & -i & 0 & 0 & 0 \\ i & 0 & -i \sqrt{\frac{3}{2}} & 0 & 0 \\ 0 & i \sqrt{\frac{3}{2}} & 0 & -i \sqrt{\frac{3}{2}} & 0 \\ 0 & 0 & i \sqrt{\frac{3}{2}} & 0 & -i \\ 0 & 0 & 0 & i & 0 \\ \end{array} \right)$$ $$S_3=\left( \begin{array}{ccccc} 2 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & -2 \\ \end{array} \right) ~~~~ $$ (with Casimir $ =6~~1\!\!1$) ?
Note the Casimir of the $\Delta(J_i)$s is not diagonal. So you must find the (sparse) Clebsch transition matrix Q such that the quadratic Casimir of the $Q\Delta(J_i) Q^{-1}$s is block diagonal, by contrast; with three diagonal blocks of size 5×5, 3×3, and 1×1, i.e. the identity matrices in 5d , 3d, and 1d, multiplied by 6, 2, and 0, respectively. (e.g., the 9,9 diagonal is 0.)