Why does reducing polynomial regression to linear regression work?
I don't really understand your problem statement. You can simply add squared terms to any linear model. Say you have $y$ and $x$ and you want to model a polynomial function, you can write a model like:

$y=\beta_0 + \beta_1 x_1 + \beta_2 x_1^2 + u$.

In matrix form this would look like 

$y=X\beta + u$.

An example with some numbers would write:

$\left( \begin{array}{rrrr}
5 \\
8 \\
7 \\
6 \\
\end{array}\right) $ = 
$\left( \begin{array}{rrrr}
1 &amp; 8&amp; 8^2  \\
1 &amp; 4 &amp; 4^2  \\
1 &amp; 3 &amp; 3^2  \\
1 &amp; 6 &amp; 6^2  \\
\end{array}\right) \beta + u.$

You can solve this like $\hat{\beta} = (X'X)^{-1} X'y$, where $\hat{\beta}$ are the coefficients of the linear regression model.