Start with $F(Z,C) = Z^2+C$, perturb this using $$f(Z,z,c)=F(Z+z,C+c)-F(Z,C)=2Zz+z^2+c$$ Claim 1 The bilinear approximation to $f$ $$g(z,c)=Az+Bc$$ with $A=2Z$ and $B=1$ has negligible relative error when $$|z| < \epsilon 2|Z| =: r_g$$ Consider chaining approximations $g_x$, $g_y$ like $$g_{y \circ x}(z,c)=A_y(A_x z + B_x c) + B_y c=(A_x A_y)z + (A_yB_x+B_y)c$$ Claim 2 this has negligible relative error compared to the original composition of $f$ functions when $$|z|< \min\left(r_x, \max\left(0,\frac{r_y-|B_x| |c|}{|A_x|}\right)\right)=:r_{y \circ x}$$
However $$f(Z_y, f(Z_x, z,c),c)=2Z_y(2Z_xz +z^2+c) +(2Z_xz+z^2+c)^2+c\\=(4Z_yZ_x)z + (4Z_x^2 +2Z_y) z^2+(2Z_y+1)c + \text{higher order terms}$$ when bilinearly approximated has negligible relative error when $$|z|<\epsilon \frac{4 |Z_y||Z_x|}{4|Z_x|^2+2|Z_y|}=:r_f$$ but (assuming $c=0$) $$r_f = \epsilon \frac{2|Z_x|}{1+4|Z_x|^2/(2|Z_y|)} < \epsilon 2|Z_x| = r_x$$ and $$r_f = \epsilon \frac{2|Z_y|}{2|Z_x| + 2|Z_y|/(2|Z_x|)} < \frac{r_y}{2|Z_x|} $$ so the required $r_f$ is less than the estimated $r_{y \circ x}$, which is a problem: i.e., claim 2 is false, for two steps, and it likely gets worse when repeating composition many times.
I have heard that Taylor's theorem can be used for rigourous error bounds. Would it be possible to get bounds for the composition of bilinear $g$ approximating $f$, by considering the composition of quadratic forms $h$ approximating $f$? Or does the full Taylor series for compositions of $f$ need to be carried along? What would be $r_{y\circ x}$ for $h(z,c)=Az +Bc +Dz^2+Ezc+Fc^2$ given $h_x,h_y,r_x,r_y$?
Actually, gerrit on fractalforums.org convinced me that claim 2 is good enough when using finite precision (e.g. floating point on a computer):
$|z| < r$ means that the quadratic term is less than rounding error in the non-linear function, so the rounding error in the output is larger than any error from non-linearity.
If $|z_x| < \min(r_x, (r_y - |B_x| |c|) / |A_x|)$, then $|z_y| \le |A_x| |z_x| + |B_x| |c| < r_y$, and $|z_x| < r_x$, so all the non-linearities are less than rounding errors at each stage.
Rounding errors can be rigourously handled via interval arithmetic.