Below we include some sample point data in a DataTable component and associate it with the variable $\mathrm{data}$.

Alternatively, this data could be imported from an external source using the Import command.

The general form of a bivariate polynomial of total degree $n$ is given by:

$f:=n\to \sum _{j\=0}^n\left(\sum _{i\=0}^n{a}_{i\,j}{x}^i{y}^j\right)$

For example, the general form for a bivariate quadratic is:

$f\left(2\right)$

Next, we choose the order of the bivariate polynomial which we will fit to the points. Increasing this value will refine the fit.

$n ≔ 3$

$\mathrm{poly} ≔ \mathrm{unapply}\left(f\left(n\right)\,x\,y\right)$

Separate and normalize the data

$\mathrm{with}\left(\mathrm{LinearAlgebra}\right)\:$

$X ≔ \mathrm{Normalize}\left(\mathrm{data}\left[..\, 1\right]\, \mathrm{Euclidean}\right)$

$Y ≔ \mathrm{Normalize}\left(\mathrm{data}\left[..\, 2\right]\, \mathrm{Euclidean}\right)$

$Z ≔ \mathrm{data}\left[..\, 3\right]$

$\mathrm{nRows} ≔ \mathrm{RowDimension}\left(\mathrm{data}\right)$

Define objective function:

$\mathrm{sse} ≔ \mathrm{add}\left({\left(\mathrm{poly}\left(X\left[i\right]\, Y\left[i\right]\right)-Z\left[i\right]\right)}^2\, i \= 1 .. \mathrm{nRows}\right) \:$

Minimize the objective function:

$\mathrm{results} ≔ \mathrm{Optimization}:-\mathrm{Minimize}\left(\mathrm{sse}\, \mathrm{iterationlimit} \= 10000000\right)$

Assign the values corresponding to the minimum value to the parameters:

$\mathrm{assign}\left(\mathrm{results}\left[2\right]\right)$

Original Data:

$\mathrm{p1} ≔ \mathrm{plots}:-\mathrm{pointplot3d}\left(\left[X\,Y\,Z\right]\,\mathrm{color}\=\mathrm{black}\,\mathrm{symbol}\=\mathrm{solidsphere}\,\mathrm{symbolsize}\=10\right)$

Best Fit Surface:

$\mathrm{p2}:=\mathrm{plot3d}\left(\mathrm{poly}\left(x\,y\right)\,x\=\min \left(X\right)..\max \left(X\right)\,y\=\min \left(Y\right)..\max \left(Y\right)\,\mathrm{style}\=\mathrm{patchnogrid}\,\mathrm{transparency}\=0.5\right)$

$\mathrm{plots}:-\mathrm{display}\left(\mathrm{p1}\,\mathrm{p2}\right)$