Multivariate Quadratic Regression I would like to make a polynomial regression, but for multivariate input data. In the univariate case, one can write polynomial regression as a multivariate linear regression problem and can come up with the closed form for ordinary least squares of
$$
\begin{pmatrix}a\\b\\c\end{pmatrix} = (\mathbf X^T \mathbf X)^{-1} \mathbf X^T \mathbf Y
$$
(see e.g. Equations For Quadratic Regression or https://en.wikipedia.org/wiki/Polynomial_regression).
However, in my case, the quadratic regression is multivariate, so
$$
\min_{a,b,C} \sum_{i=1}^N y_i - (a + b^T\cdot x_i + x_i^T\cdot C\cdot x_i)^2
$$
where $C$ is a symmetric matrix, $b$ and $x_i$ are vectors, $y_i$ and $a$ are scalars and $N$ is the number of samples (we can assume we have enough samples to have an overdetermined system).
Does a closed form exist here as well, and if so, what does it look like?
If not, how do I do the regression? Obviously I could use regular optimization methods, like BFGS, with the constraints, that C is symmetric, but that is not as efficient as I would hope for.
 A: 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}
$$
A: 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.
