I would like to fit 2D data (a set of $x_i$, $y_i$, $z_i = f(x_i, y_i)$) where my function is a product of two 1D polynomials of $n$th order ($n>1$), so it is factorizable.
$$f(x,y) = (a_0 + a_1x + a_2x^2 + ... + a_nx^n)(1 + b_1y + b_2y^2 + ... + b_nx^n)$$
I know that I can solve this using any optimization routine that uses non-linear least squares. Is there a way to use a linear least square method to fit such a factorizable polynomial?
One can expand the polynomial and formulate the problem as
$$ \begin{pmatrix} 1 & x_1 & y_1 & x_1y_1 & x_1^2 & y_1^2 & x_1^2y_1 & x_1y_1^2 & \cdots \\ &&&& \vdots &&&& \end{pmatrix} \mathbf{X} = \begin{pmatrix} f(x_1, y_1) \\ \vdots \end{pmatrix} $$ and solve for $\mathbf{X}= (a_0, a_1, a_0b_1, a_1b_1, \cdots)^T$.
Then I have the problem, that the system of equations that define $a_i$ and $b_i$ is overdetermined ($\mathbf{X}$ has length $(n+1)^2$, but I have only $2n+1$ independent variables).
Constrained least squares (using e.g. $\mathbf{L}\mathbf{X} = \mathbf{d}$) also does not seem to work, since the constraints are non-linear.
Is there any (tricky?) way to use linear least squares for this type of problem?
Example of such data that I try to fit to ($n=4$): $x = [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , -13.5 , -10.125 , -6.75 , -3.375 , 0. , 3.375 , 6.75 , 10.125 , 13.5 , 13.510486 , 13.50262 , 13.50262 , 13.510486 , 6.755392 , 6.7513475, 6.7513475, 6.755392 , -6.755392 , -6.7513475, -6.7513475, -6.755392 , -13.510486 , -13.50262 , -13.50262 , -13.510486 ]$ $y = [ 0. , 0. , 2.3 , 1.15, 0. , -1.15, -2.3 , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 2.3 , 1.15, -1.15, -2.3 , -2.3 , -1.15, 1.15, 2.3 , 2.3 , 1.15, -1.15, -2.3 , -2.3 , -1.15, 1.15, 2.3 ]$ $f(x, y) = [1.009518 , 0.99437314, 1.0034443 , 1.0105053 , 1.0142992 , 1.0154775 , 0.9043408 , 0.9617672 , 1.000238 , 1.0179688 , 1.0113076 , 0.98180985, 0.9331824 , 0.8696933 , 0.79686165, 0.7879401 , 0.7935547 , 0.7991753 , 0.79890484, 0.93638456, 0.935677 , 0.92781633, 0.9205955 , 0.9865708 , 0.99512374, 1.0054537 , 1.0065441 , 0.9099518 , 0.90881765, 0.9004598 , 0.892026 , 1.013272 ]$
The method below is not verry accurate but might be sufficient. For example the fitting is not bad with your data as shown below.
$$f(x,y) = g(x)h(y)\quad \begin{cases} g(x)=(a_0 + a_1x + a_2x^2 + ... + a_nx^n) \\ h(y)=(1 + b_1y + b_2y^2 + ... + b_ny^n) \end{cases}$$ We approximate $g(x)$ with the function : $$g(x)\simeq \frac{1}{(c_0 + c_1x + c_2x^2 + ... + c_mx^m)}$$ Note that $m$ can be different from $n$. $$f(x,y)(c_0 + c_1x + c_2x^2 + ... + c_mx^m) \simeq (1 + b_1y + b_2y^2 + ... + b_nx^n)$$ $$c_0f + c_1xf + c_2x^2f + ... + c_mx^mf-b_1y - b_2y^2 - ... - b_ny^n\simeq 1$$ This is a linear equation wrt the coefficients. Thus $c_0,c_1,...,c_m,b_1,b_2,b_n$ can be approximately computed thanks to linear regression.
Now knowing the approximate values of $b_1,b_2,...,b_n$ we can compute
$$H_k=b_1y_k + b_2y_k^2 + ... + b_ny_k^n \quad \text{from }k=1 \text{ to } N=\text{number of points.}$$
The original equation becomes : $$f_k\simeq a_0H_k + a_1H_kx_k + a_2H_kx_k^2 + ... + a_nH_kx_k^n$$ This is a linear equation wrt the coefficients. Thus $a_0,a_1,...,a_n$ can be approximately computed thanks to linear regression.
So, approximate values of $a_0,a_1,...,a_n,b_1,...,b_n$ are obtained.
Numerical example with your data :
The calculus was made with $m=4.$
$$c_0f + c_1xf + c_2x^2f + c_3x^3f + c_4x^4f-b_1y - b_2y^2 - b_3y^3 - b_4y^4\simeq 1$$ The linear regression gives :
$c_0=1.027432 \: ; \: c_1=0.004248 \:;\: c_2=-1.322825 \:10^{-4} \:;\: c_3=9.926459 \:10^{-6} \:;\: c_4=4.165036 \:10^{-6}$
$b_1=-0.005547\:;\:b_2=-0.010356\:;\:b_3=0.002587\:;\:b_4=0.002165$
Then the $H_k=b_1y_k + b_2y_k^2 + ... + b_ny_k^n$ were computed.
The linear regression from $\quad f_k\simeq a_0H_k + a_1H_kx_k + a_2H_kx_k^2 + a_3H_kx_k^3 + a_4H_kx_k^4\quad$ leads to : $$a_0=0.996376 \:;\: a_1=-0.006457\:;\: a_2=-0.001046\:;\: a_3=1.129896\:10^{-5}\:;\: a_4=2.196121\:10^{-6}$$
To check the result, one compute $$fcal_k=(a_0 + a_1x_k + a_2x_k^2 + a_3x_k^3 + a_4x_k^4)(1 + b_1y_k + b_2y_k^2 + b_3y_k^3 + b_4y_k^4)$$ $$RMSE=\sqrt{\frac{1}{N}\sum_{k=1}^N(f_k-fcal_k)^2}=0.0409$$ And visually :
Probably the fitting can be slightly improved with higher $m$, for example 5 or 6 , ... , but I didn't try it.
Of course non-linear regression gives a better accuracy of fitting. But this is not easy with such high number of parameters. The softwares for non-linear fitting require to give guessed initial values or guessed ranges for the parameters. In case of 9 parameters this is a big difficulty in practice. Often the calculus goes out of range or takes too long time.
This difficulty can be overcome thanks to the above method. The above method gives appoximate values of the parameters which can be used as good initial values for starting the iterative calculus of non-linear regression.