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 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 at 19:55
• Could it be that you missed two $\sum$-symbols in your design matrix? – Make42 Mar 20 at 23:00
• Could be... corrected now – Mark Fischler Mar 21 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 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 at 7:35
• @MarkFischler: Btw: I would be happy if you would take this answer and edit it into your own (basically you could copy-paste). I will accept your post as the answer: I would like to have a complete picture in the accepted answer (the mods, did not let me edit your post). – Make42 Mar 25 at 11:12
• @MarkFischler: If you explain why the second approach is so much slower, I would appreciate this too: The sums need to be done either why, so I am surprised. The only way I see an advantage of approach 1 is if I was able to only calculate the upper triangular matrix and then copy it down to the lower triangle. This would be ~2 times faster. – Make42 Mar 25 at 11:14