Gauge invariance doesn't actually force usual choice of covariant derivative? As we all know, a gauge invariant theory is of the form 
$$ \mathcal{L} = \bar{\psi} \gamma^\mu \left(  i\partial_\mu + A_\mu^a T^a\right) \psi.$$
The multiplet $\psi$ and gauge field $A_\mu = A^a_\mu T^a$ transform as follows under a gauge transformation:
$$\psi \rightarrow G \psi, \quad \quad A_\mu \rightarrow G A_\mu G^{-1} -  (\partial_\mu G)  G^{-1},$$
where $G(x)$ is an element of the gauge group.
I know that $T^a$ are supposed to be the generators of the gauge group, i.e. they are a basis for the associated Lie Algebra. However, it seems to me that actually this fact is not essential for $\mathcal{L}$ to be gauge invariant! I mean, any old matrices $T^a$ will do; we just say that under a 'gauge transformation' $A_\mu$ transforms as above. 
So I am confused. Suppose I was ignorant and all I want to do is construct a gauge invariant theory. It would seem that I should be able to take any matrix (of the right dimension) for $T^a$, or am I wrong?
If I am right then, from a constructionist point of view, what is the prime reason for choosing $T^a$ as we normally do? 
 A: Suppose that the lagrangian
$$
\newcommand{\cL}{{\cal L}}
\newcommand{\opsi}{{\overline \psi}}
\newcommand{\pl}{\partial}
 \cL=\opsi\gamma^\mu(i\pl_\mu+A_\mu)\psi
\hskip2cm
 A_\mu := \sum_a A^a_\mu T^a
\tag{1}
$$
is invariant under gauge transformations$^\dagger$
\begin{gather}
 \psi\to G\psi
\tag{2}\\
 (i\pl_\mu+A_\mu)\to G(i\pl_\mu+A_\mu)G^{-1}
\tag{3}
\end{gather}
for all $G$ in some matrix group, where the $T^a$ are matrices of the same size.  In order for this to make sense, the right-hand side must end up with the same matrices $T^a$ as the left-hand side, because the fields $A_\mu^a$ are the only things being transformed in the last equation. (The components of the matrices are just fixed coefficients in the Lagrangian, like the coefficient $m$ in a mass term.) This gives the requirements
\begin{gather}
 G\pl_\mu G^{-1} = \text{linear combination of }T^a\text{s}
\tag{4}\\
 GT^a G^{-1} = \text{linear combination of }T^a\text{s}.
\tag{5}
\end{gather}
By taking $G$ to be infinitesimally close to the identity, equation (4) implies that $G$ is generated by the $T^a$s, and equation (5) implies that the commutator of two $T^a$s must be a linear combination of $T^a$s.

$^\dagger$ The second equation in (3) expresses how $A_\mu$ transforms. The partial derivatives on the right-hand side act on both $G^{-1}$ and whatever stands to the right of $G^{-1}$, just like the partial derivative on the left-hand side acts on whatever stands to the right of the closing parenthesis.
A: Since a gauge field is spacetime-dependent, a translation in spacetime is necessarily accompanied by some change in gauge. Thus, you can think of the gauge covariant derivative ($D_\mu$) as an infinitesimal spacetime translation ($\partial_\mu$) along with an infinitesimal transformation in gauge space (any other terms appearing in $D_\mu$). Then it should hopefully be more clear that the other terms in $D_\mu$ must somehow relate to the gauge group, and cannot be written in terms of completely arbitrary matrices.
To be more concrete, recall/note that an element ${g}$ of a Lie group $G$ is given by
\begin{equation}
{g}=\exp(iA_\mu^aT^a)\,,
\end{equation}
where $A_\mu^a$ are the continuous parameters of $G$, and ${T}^a$ its generators. We can use a Taylor expansion to write this as:
\begin{equation}
{g}={I}+\sum_{n=1}^\infty \frac{1}{n!}\,(iA_\mu^a{T}^a)^n={I}+iA_\mu^a{T}^a+\mathcal{O}\left((A_\mu^a)^2\right),
\end{equation}
where $I$ is the identity element. For the infinitesimal transformation one takes the leading term in the expansion ($iA_\mu^a{T}^a$), and this is what you see in your example:
\begin{equation}
\mathcal{L}={\bar\Psi}i\gamma^\mu {D}_\mu {\Psi}={\bar\Psi}i\gamma^\mu (\partial_\mu+iA_\mu^a {T}^a) {\Psi}\,.
\end{equation}
