Gauge theory formalism Source:
Chapter 11 of Freedman and Van Proeyen’s supergravity textbook

An infinitesimal symmetry transformation is determined by 
1) a parameter, call it $\epsilon^A$, and 
2) an operation, call it $\delta_\epsilon$.
The operation $\delta_\epsilon$
1) depends linearly on the parameter $\epsilon^A$, and 
2) acts on fields, i.e. $\delta_\epsilon(\phi^i)$.
For some global symmetry, $\epsilon^A$ does not depend on the spacetime $x^\mu$.
Another way to say "$\delta_\epsilon$ depends linearly on the parameter $\epsilon^A$,'' is to write 
$$\delta_\epsilon =\epsilon^A T_A$$
where the $T_A$ are some operations on fields.
Let $\{(t_A)^i{}_j\}$ be the matrix generators of a representation of some Lie algebra.
This Lie algebra (LA) is defined by $[t_A,t_B]=f_{AB}{}^C t_C$.
The action of $T_A$ on the fields is defined with the LA basis elements,
$$ T_A(\phi^i)=-(t_A)^i{}_j \phi^i $$
So then we have 
\begin{eqnarray*}
\delta_\epsilon(\phi^i) &=& \epsilon^A T_A(\phi^i) \\
&=& -\epsilon^A (t_A)^i{}_j (\phi^j)
\end{eqnarray*}
Then the product of two symmetry transformations reads,
\begin{eqnarray*}
\delta_{\epsilon_1}\delta_{\epsilon_2}(\phi^i) &=& \epsilon_1{}^A T_A(\epsilon_2{}^B T_B\phi^i) \\
&=& \epsilon_1{}^A T_A(-\epsilon_2{}^B (t_B)^i{}_j \phi^j) \\
&=& -\epsilon_1{}^A \epsilon_2{}^B (t_B)^i{}_j T_A \phi^j \\
&=& -\epsilon_1{}^A \epsilon_2{}^B (t_B)^i{}_j (-(t_A)^j{}_k \phi^k) \\
&=& \epsilon_1{}^A \epsilon_2{}^B (t_B)^i{}_j (t_A)^j{}_k \phi^k \\
\end{eqnarray*}

The authors then go on to state the commutator, which I am concerned with,

Notably,
$\epsilon_1{}^A$ and  $\epsilon_2{}^B$ are numbers, and are so are commutative, and
$(t_B)^i{}_j,$ $(t_A)^j{}_k,$ and  $\phi^k$ are matrices, and so are associative.
Thus, how is $[\delta_{\epsilon_1},\delta_{\epsilon_2}](\phi^i) = \delta_{\epsilon_1}\delta_{\epsilon_2}(\phi^i) - \delta_{\epsilon_2}\delta_{\epsilon_1}(\phi^i)$
not equal to zero?
As far as I can tell, $\delta_{\epsilon_1}\delta_{\epsilon_2}(\phi^i)=\delta_{\epsilon_2}\delta_{\epsilon_1}(\phi^i).$

Disclaimer:
That formula for the product of symmetry transformations on the fields is my work. The authors have, 

which seems the same as mine but I just want to be cautious.
Cheers
 A: I think its mostly a notational confusion. You derivation for the $\delta$ commutator does lead to the same conclusion as the authors.
$$
\begin{eqnarray*}
\delta_1\delta_2
&=&\epsilon_1{}^A \epsilon_2{}^B (t_B)^i{}_j (t_A)^j{}_k \phi^k \\
\delta_2\delta_1
&=&\epsilon_2{}^A \epsilon_1{}^B (t_B)^i{}_j (t_A)^j{}_k \phi^k \\
&=& \epsilon_1{}^B \epsilon_2{}^A (t_B)^i{}_j (t_A)^j{}_k \phi^k \\
&=& \epsilon_1{}^A \epsilon_2{}^B (t_A)^i{}_j (t_B)^j{}_k \phi^k 
\end{eqnarray*}
$$
where the last step follows from $A,B$ being dummy vars.
$$
\begin{eqnarray*}
[\delta_1,\delta_2]
&=&  \epsilon_1{}^A \epsilon_2{}^B\left( (t_B)^i{}_j (t_A)^j{}_k-(t_A)^i{}_j (t_B)^j{}_k\right) \phi^k   \\
&=& \epsilon_1{}^A \epsilon_2{}^B\left( (t_B t_A)^i{}_k-(t_At_B)^i{}_k\right) \phi^k \\
&=& \epsilon_1{}^A \epsilon_2{}^B([t_B,t_A] )^i{}_k \phi^k\\
&=& \epsilon_1{}^A \epsilon_2{}^B(-f_{AB}{}^Ct_C )^i{}_k \phi^k\\
&=& \epsilon_1{}^A \epsilon_2{}^Bf_{AB}{}^C(-t_C )^i{}_k \phi^k\\
&=& \epsilon_1{}^A \epsilon_2{}^Bf_{AB}{}^CT_C \phi^i\\
\end{eqnarray*}
$$
A: Thank you to @lineage for the prompting!
\begin{eqnarray*}
\delta_{\epsilon_1}\delta_{\epsilon_2}(\phi^i) &=& \epsilon_1{}^A \epsilon_2{}^B (t_B)^i{}_j (t_A)^j{}_k \phi^k \\
&\neq& \epsilon_1{}^A \epsilon_2{}^B (t_A)^i{}_j (t_B)^j{}_k \phi^k
\end{eqnarray*}
since matrix multiplication is not commutative.
Thus the two terms in the commutator are not identical.
