Suppose we have an action $S=S(a,b,c)$ which is a functional of the fields $a,\, b,\,$ and $c$. We denote the variation of $S$ wrt to a given field, say $a$, i.e. $\frac{\delta S}{\delta a}$, by $E_a$.

Then $S$ is gauge invariant when

$$\delta S = \delta a E_a + \delta b E_b +\delta c E_c = 0 \tag{1}$$

This gives

$$\delta c = - (\delta a E_a + \delta b E_b)/E_c \tag{2}$$

From the above well-defined equation $\delta c$ can be obtained for arbitrary $\delta a$ and $\delta b$. It is not necessary to have a relation b/w $\,\delta a$ and $\delta b$. So there is an infinite number of gauge-transformations here. These gauge transformations vanish when the equations of motion hold, for e.g. $\delta c$ in $(2)$ vanishes when $E_a=0,\,\&\,E_b=0$. In other words, these gauge transformations vanish on-shell.

As per sec. 3.1.5 of Ref. 1, such gauge transformations are called "trivial gauge transformations". According to theorem 3.1 such gauge transformations can always be expressed as $\delta y^i=\varepsilon^{ij}\frac{\delta S}{\delta y^j}$ where the action $S$ is a functional of the fields $y^i$ and $\varepsilon^{ij}=-\varepsilon^{ji}$.

I am having trouble obtaining the components of the matrix $\varepsilon^{ij}$ for the above example. If $y^1=a,\,y^2=b,$ and $y^3=c\,$ then from $\delta y^1=-(\delta bE_b+\delta cE_c)/E_a\,$ I read $\varepsilon^{12}=(-\delta b)/E_a$. From $\delta y^2=-(\delta aE_a+\delta cE_c)/E_b\,$ I read $\varepsilon^{21}=(-\delta a)/E_b$. But these two components are not simply (-1) times each other as they should be. What's going wrong?

References:

1. M. Henneaux & C. Teitelboim, Quantization of Gauge Systems, 1994; subsection 3.1.5.