This question shows how to least square fit a function to some data points: Creating a function incrementally
That question/answer show how to fit a function of the form $y=f(x)$ using $(x,y)$ pairs of data points.
That can be generalized to surfaces and hyper-volumes by instead of having powers of $x$, having the permutations of all the powers of each of the arguments. For example, the function $z=f(x,y)=Axy+Bx+Cy+D$ would have an $A^TA$ matrix and $A^Tz$ vector that consisted of the combinations of the powers of $x$ and $y$: $x^1y^1,x^1y^0,x^0y^1,x^0y^0$ otherwise written as $xy,x,y,1$.
The pattern continues for higher dimension functions such as $w=f(x,y,z)$ and higher.
That all makes well enough sense to me but now I'm wondering, what if I have a vector function such as $(z,w)=f(x,y)$?
Is there a way, given $(x,y)$ points to do a least square fit function to those points?
If there isn't, are there any other decent / common methods to fit a function of that form to points?
Consider the function $f\colon\mathbb{R}^{\nu}\mapsto\mathbb{R}^{\mu}$. We have a sequence of $m$ measurements $$ \left\{ \left( x_{1_{k}}, x_{2_{k}}, x_{3_{k}}, \dots, x_{\nu_{k}} \right), \left( y_{1_{k}}, y_{2_{k}}, y_{3_{k}}, \dots, y_{\mu_{k}} \right) \right\}_{k=1}^{m} $$ which we choose to approximate with a fit of order $d$: $$ f \left( x_{1}, x_{2}, x_{3}, \dots, x_{\nu} \right) = c_{000\dots0} + c_{100\dots0}x_{1} + c_{010\dots0}x_{2} + c_{001\dots0}x_{3} + c_{000\dots1}x_{\nu} + c_{200\dots0}x_{1}^{2} + c_{110\dots0}x_{1}x_{2} + c_{101\dots0}x_{1}x_{3} + c_{000\dots2}x_{\nu}^{2} + \dots + c_{000\dots d} x_{\nu}^d $$ The linear system is $$ \begin{align} \mathbf{A} c &= F \\ \left[ \begin{array}{cccc} 1 & x_{1_{1}} & x_{2_{1}} & x_{3_{1}} & \dots & x_{\nu_{1}}^{d} \\ 1 & x_{1_{2}} & x_{2_{2}} & x_{3_{2}} & \dots & x_{\nu_{2}}^{d} \\ \vdots & \vdots & \vdots & \vdots & & \vdots\\ 1 & x_{1_{m}} & x_{2_{m}} & x_{3_{m}} & \dots & x_{\nu_{m}}^{d} \end{array} \right] % \left[ \begin{array}{cccc} c_{000\dots0_{1}} & c_{000\dots0_{2}} & \dots & c_{000\dots0_{\mu}} \\ c_{100\dots0_{1}} & c_{100\dots0_{2}} & \dots & c_{100\dots0_{\mu}} \\ c_{010\dots0_{1}} & c_{010\dots0_{2}} & \dots & c_{010\dots0_{\mu}} \\ c_{001\dots0_{1}} & c_{001\dots0_{2}} & \dots & c_{001\dots0_{\mu}} \\ \vdots & \vdots && \vdots \\ c_{000\dots1_{1}} & c_{000\dots1_{2}} & \dots & c_{000\dots1_{\mu}} \\ c_{200\dots0_{1}} & c_{200\dots0_{2}} & \dots & c_{200\dots0_{\mu}} \\ \vdots & \vdots && \vdots \\ c_{000\dots d_{1}} & c_{000\dots d_{2}} & \dots & c_{000\dots d_{\mu}} \end{array} \right] &= \left[ \begin{array}{c} % f \left( \bar{x}_{1} \right)_{1} & f \left( \bar{x}_{1} \right)_{2} & \dots & f \left( \bar{x}_{1} \right)_{\mu} \\ % f \left( \bar{x}_{2} \right)_{1} & f \left( \bar{x}_{2} \right)_{2} & \dots & f \left( \bar{x}_{2} \right)_{\mu} \\ % \vdots & \vdots && \vdots \\ % f \left( \bar{x}_{\nu} \right)_{1} & f \left( \bar{x}_{\nu} \right)_{2} & \dots & f \left( \bar{x}_{\nu} \right)_{\mu} \\ % \end{array} \right] \\ \end{align} % $$
Notation: $\bar{x}_{2}$ represents all $x$ values in the second measurement; $f\left(\bar{x}_{2}\right)_{3}$ is the third component of the vector function $f$ evaluated at the second measurement location. $$ \mathbf{A}\in\mathbb{C}^{m\times n}, \quad c\in\mathbb{C}^{n \times \mu}, \quad F\in\mathbb{C}^{m \times \mu}. $$