How to do regression using polynomials of the form $a_n x^n + a_{n-1} x^{n-1} y + \cdots+ a_0 y^n+ b_{n-1}x^{n-1}+b_{n-2}x^{n-2}y+\cdots=c$, given some data points $(x_i,y_i)$ (the number of data points is much greater than $2n+2$)?
The sources on the internet only seem to suggest the method for polynomial functions of the form $y=f(x)$.
On
Just do a linear regression like you normally would with any basis of functions: calculate the $x^ky^l$, put in the columns of a matrix $\bf \Phi$ in some order, then pose
$${\bf v_o}= \min_{\bf v}\|{\bf \Phi v}-c{\bf 1}\|_2^2$$
If you allow constant term ($k=l=0$), there will exist a trivial solution which you will need to dodge somehow.
For example by regularization on the coefficients of the polynomial.
The value at position $m$ in the vector $\bf v$ will be the coefficient for whichever of the $x^ky^l$ you stuffed as $m$th column in the $\bf \Phi$.
Assuming that you have $p$ data points $(x_i,y_i)$ and that what is measured is $y$, in the least square sense, I do no thing that this is feasible.
One think you could do is to consider that you face a surabondant set of equations. So, just define $z_i=c \,\,\forall i$ and use multilinear regression or matrix calculations to get the $a$'s and $b$'s.