Multivariate Quadratic Regression
Solution 1:
You can do multi-variate quadratic regression in the usual way. Let's label the row (and column) indices of the design matrix $A$, and the row index of the value vector $b$, by index $s(\{p_1, p_2, p_3, \cdots\})$ which pertains to the coefficient of $x_i^{p_1}x_2^{p_2}\cdots$. For example, the row labeled $s(\{ 1, 0, 2\})$ will be the row pertaining to the coefficient of $x_1x_3^2$.
Then the elements of $A$ are calculated as $$ A_{s(\{p_1, p_2, p_3, \cdots\}),s(\{q_1, q_2, q_3, \cdots\})} = \sum x_1^{p_1+q_1} x_2^{p_2+q_2} x_3^{p_3+q_3} \cdots $$ and the elements of $b$ are $$ b_{s(\{p_1, p_2, p_3, \cdots\})} = \sum y\,x_1^{p_1} x_2^{p_2} x_3^{p_3} \cdots $$ where of course all the sums are taken over the set of data points.
For example, for a 2-variable quadratic fit $y = a + bu + cv + du^2 + e uv + fv^2$ you need to solve $$ \pmatrix{N &\sum u_i &\sum v_i & \sum u_i^2 & \sum u_iv_i & \sum v_i^2 \\ \sum u_i & \sum u_i^2 & \sum u_i v_i & \sum u_i^3 & \sum u_i^2v_i & \sum u_i v_i^2 \\ \sum v_i & \sum u_iv_i & \sum v_i^2 & \sum u_i^2v_i & \sum u_iv_i^2 & \sum v_i^3 \\ \sum u_i^2 & \sum u_i^3 & \sum u_i^2 v_i & \sum u_i^4 & \sum u_i^3v_i & \sum u_i^2 v_i^2 \\ \sum u_iv_i & \sum u_i^2v_i & \sum u_i v_i^2 & \sum u_i^3v_i & \sum u_i^2v_i^2 & \sum u_i v_i^3 \\ \sum v_i^2 & \sum u_iv_i^2 & \sum v_i^3 & \sum u_i^2v_i^2 & \sum u_iv_i^3 & \sum v_i^4 } \pmatrix{a\\b\\c\\d\\e\\f} =\pmatrix{\sum y_i \\ \sum y_i u_i \\ \sum y_iv_i \\ \sum y_iu_i^2\\ \sum y_iu_iv_i \\ \sum y_iv_i^2} $$
Solution 2:
Disclaimer: Approach 1 is from Mark Fischler, but I want to reference the approach in my second approach and I need the labels under the matrices for referencing, so I restate the approach. Apparently, adding the second approach to Mark's answer is not wanted by the moderators.
Approach 1
You can do multi-variate quadratic regression in the usual way. Let's label the row (and column) indices of the design matrix $A$, and the row index of the value vector $b$, by index $s(\{p_1, p_2, p_3, \cdots\})$ which pertains to the coefficient of $x_i^{p_1}x_2^{p_2}\cdots$. For example, the row labeled $s(\{ 1, 0, 2\})$ will be the row pertaining to the coefficient of $x_1x_3^2$.
Then the elements of $A$ are calculated as $$ A_{s(\{p_1, p_2, p_3, \cdots\}),s(\{q_1, q_2, q_3, \cdots\})} = \sum x_1^{p_1+q_1} x_2^{p_2+q_2} x_3^{p_3+q_3} \cdots $$ and the elements of $b$ are $$ b_{s(\{p_1, p_2, p_3, \cdots\})} = \sum y\,x_1^{p_1} x_2^{p_2} x_3^{p_3} \cdots $$ where of course all the sums are taken over the set of data points.
For example, for a 2-variable quadratic fit $y = a + bu + cv + du^2 + e uv + fv^2$ you need to solve $$ \underbrace{\pmatrix{N &\sum u_i &\sum v_i & \sum u_i^2 & \sum u_iv_i & \sum v_i^2 \\ \sum v_i & \sum u_iv_i & \sum v_i^2 & \sum u_i^2v_i & \sum u_iv_i^2 & \sum v_i^3 \\ \sum u_i^2 & \sum u_i^3 & \sum u_i^2 v_i & \sum u_i^4 & \sum u_i^3v_i & \sum u_i^2 v_i^2 \\ \sum u_iv_i & \sum u_i^2v_i & \sum u_i v_i^2 & \sum u_i^3v_i & \sum u_i^2v_i^2 & \sum u_i v_i^3 \\ \sum v_i^2 & \sum u_iv_i^2 & \sum v_i^3 & \sum u_i^2v_i^2 & \sum u_iv_i^3 & \sum v_i^4 }}_{\mathbf A} \pmatrix{a^*\\b^*\\c^*\\d^*\\e^*\\f^*} = \underbrace{ \pmatrix{\sum y_i \\ \sum y_i u_i \\ \sum y_iv_i \\ \sum y_iu_i^2\\ \sum y_iu_iv_i \\ \sum y_iv_i^2} }_{\mathbf b} $$
where $a^*, b^*, c^*, d^*, e^*, f^*$ are the optimal values of $a, b, c, d, e, f$ after the quadratic fit.
Approach 2
Alternatively we can consider
\begin{align} \mathbf Y &= \mathbf X\cdot\pmatrix{a\\\dots\\f}% \\ \underbrace{\pmatrix{y_{1}\\y_{2}\\y_{3}\\\vdots \\y_{n}}}_{\mathbf Y} &= \underbrace{\pmatrix{ 1&u_1&v_1&u_1^2 & u_1v_1 & v_1^2\\ 1&u_2&v_2&u_2^2 & u_2v_2 & v_2^2\\ 1&u_3&v_3&u_3^2 & u_3v_3 & v_3^2\\ \vdots &\vdots &\vdots &\vdots &\vdots&\vdots \\ 1&u_n&v_n&u_n^2 & u_nv_n & v_n^2\\ }}_{\mathbf X} \cdot \pmatrix{a\\b\\c\\d\\e\\f} \end{align}
We can use this to use the regular formula from https://en.wikipedia.org/wiki/Polynomial_regression for Ordinary Least Squares and get
\begin{align} \pmatrix{a^*\\b^*\\c^*\\d^*\\e^*\\f^*} = {(\underbrace{\mathbf X^{\mathsf T}\cdot\mathbf X}_{\mathbf A} )}^{-1} \cdot \underbrace{\mathbf{X}^{\mathsf T}\cdot \vec {y}}_{\mathbf b} \end{align}
Calculate your original quadratic function
You can simply
\begin{align} \alpha^* &= a^*\\ \mathbf \beta^* &= \pmatrix{b^*\\c^*}\\ \mathbf \Gamma^* &= \pmatrix{d^*&e^*\\e^*&f^*} \end{align}
for your original problem
$$ \min_{A,B,C} \sum_{i=1}^N y_i - (\alpha + \mathbf \beta^T\cdot x_i + x_i^T\cdot \mathbf \Gamma\cdot x_i)^2 $$
where $\alpha$ is a scalar, $\mathbf \beta$ is a vector and $\mathbf \Gamma$ is a matrix.