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

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.

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

• Is there a way to write the design matrix $A$ as the result of matrix multiplications, instead of calculating each element independently? Then the calculation of the design matrix would be much more efficient in the respective computer systems like Matlab or numpy. – Make42 Mar 20 '19 at 19:47
• So $C = \begin{pmatrix}d&e\\e&f\end{pmatrix}$, $b^T = (b,c)$ (first my vector than your scalars) and my $a$ equals your $a$? – Make42 Mar 20 '19 at 19:55
• Could it be that you missed two $\sum$-symbols in your design matrix? – Make42 Mar 20 '19 at 23:00
• Could be... corrected now – Mark Fischler Mar 21 '19 at 2:25
• The calculation of the matrix is incredibly efficient in any sort of procedural language, where you simply accumulate, as each point is processed, the various components. I've seen Matlab code that tried to use the powerful matrix syntax by artificially creating a column matrix out of each point - it is both ugly and very inefficient. – Mark Fischler Mar 21 '19 at 2:28

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}

You can simply

\begin{align} \alpha^* &= a^*\\ \mathbf \beta^* &= \pmatrix{b^*\\c^*}\\ \mathbf \Gamma^* &= \pmatrix{d^*&e^*\\e^*&f^*} \end{align}

$$\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.
• The second approach is fine mathematically, but is $6$ to $10$ times more calculation (less efficient) than the first. – Mark Fischler Mar 23 '19 at 7:35