Non-linear optimization of a family of functions

I'd like some advice about the following problem (resolution methods, problem category, etc.). The context is about the coregistration of several images all-in-once ($$n$$ images, $$f_{i,j}$$ the function that maps coordinates $$(x_i, y_i)$$ to $$(x_j, y_j)$$.

Let $$f_{i,j}^k(x, y)$$ a polynomial function of order $$k$$ from $$\mathbb{R}^2$$ to $$\mathbb{R}^2$$, with $$i, j \in [1, n]$$.

For example, if $$k=0$$, such functions can be expressed as $$f_{i,j}^0(x, y) = (a_{ij}^0, a_{ij}^1)$$. Or if $$k = 1$$: $$f_{i,j}^1(x, y) = (a_{ij}^0 + b_{ij}^0 x + c_{ij}^0 y, a_{ij}^1 + b_{ij}^1 y + c_{ij}^1 y)$$.

Let's assume $$k$$ is fixed, it will be omitted in the following.

This family of functions should respect as much as possible following the property: $$f_{i,j} = f_{i, k}(f_{k, j})$$.

I want to use this property to decrease the number of unknowns and obtain robust estimates of coefficients. Moreover, I have a large set of matches that can be expressed as such: $$v = f_{i,j}(u)$$.

My problem is the following: how can I compute robustly the functions $$f_{1,i}$$ ?

If $$k = 0$$, the problem is quite simple: $$f_{i,j} = f_{i, k}(f_{k, j}) = f_{i, k} + f_{k, j}$$ so i can write an overdetermined linear system and solve it with a Moore-Penrose inverse or an algorithm such as RANSAC.

If k = 1 or 2, I don't really known how to proceed. I think I could try to design a custom convergence scheme with a predetermined equation resolution order and some iterations to converge.

As an example, if I solve for $$f_{1,2}$$, then to get $$f_{1,3}$$ I can use my matches of the form $$v = f_{1,3}(u)$$ but also the matches such as $$v = f_{2,3}(u) \Rightarrow f_{1, 2}(v) = f_{1,3}(u)$$