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\cend{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^Tcdot x_i + x_i^Tcdot Ccdot 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.

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}
$$

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 Xcdotpmatrix{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}cdotmathbf 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^Tcdot x_i + x_i^Tcdot mathbf Gammacdot x_i)^2
$$

where $alpha$ is a scalar, $mathbf beta$ is a vector and $mathbf Gamma$ is a matrix.