Since they consider $A'$ and $\lambda'$ to be local functions, they can only be those constructed form $A,\partial A, \lambda, \partial \lambda$ and possibly higher order derivatives. Since there's only the constraint on the correspondences on gauge transformation, they probably guessed some forms of combinations of those variables, and find one possible solutions. They did not claim the map is unique, so you are probably right in that the single equation system is under-constrained to completely determine the map.
By treating $A$ and $\partial A$ as independent variables, and guessing $A'_i(A)=f_{ij}(A) A_j$ and $\lambda'(A,\lambda)=g_i(A,\lambda)A_i$. Setting $\lambda=0$ and $A=0$, the constraint becomes $f_{ij}(A)\partial_j \lambda-\partial_i (g_j(A,\lambda) A_j)=-\frac12\theta^{kl}\{\partial_k\lambda,\partial_l A_i\}.$ If we further guess $g_j(A,\lambda)A_j=g_{ij}\partial_i\lambda A_j$, replacing $\partial_j\lambda$ by $A_j$ leads to something mixture of $\theta^{kl}$ and $g_{kl}$ in the expression of $A'$ which now depends on $A$ and $\partial A$.
Some other choices on setting $\lambda$, $A$ and $\partial \lambda$, $\partial A$ can then be used to simultaneous determine $g_{kl}$ and $f_{ij}(A)$. So my guess is that they did some trial-and-error and figured out one (probably non-unique) map.
I guess it wouldn't be too bizarre to guess a solution that would depend only on the pointwise values and the first-order derivatives. And if $\partial A$, $A$ and $\partial r$ are treated as independent choices for the mapping, $A'$ would be scaling linearly with $A$ and linearly with $\partial A$, whereas $r'$ scales linearly with $A$ and $\partial r$. What's left is then to figure out the undetermined coefficients that must be in the order of $O(\theta)$.