Least square fitting multiple values 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? 
 A: 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}.
$$
