In Scikit-learn I can use polynomial features to create polynomial linear regression models. Scikit-learn transforms my original data as follows. If I had a datapoint [2], it would return [2,4,8] with 3 basis functions $(x, x^2, x^3).$ This means that it is increasing the dimensionality of the model, indeed because these monomials form a basis for the vector space of cubics. So the resulting model is a linear combination of monomials.

My question therefore is how do basis functions relate to the dimensionality of the model? Precisely which vector space do they form a basis for? How do we construct the linear regression model taking into account that vector space?


1 Answer 1


The model space

As always in a linear model, with a $n\times k$ model matrix $X$ for $n$ observations and $k$ features, the space is the set of all possible vectors $X\beta$ where $\beta$ is the collection of $k$ coefficients (the "parameters").

Geometrically, the parameter $\beta\in\mathbb{R}^k$ is a vector and $X$ represents a linear transformation from $\mathbb{R}^k$ to $\mathbb{R}^n.$ The model space is the image of $X.$ It's a linear subspace of $\mathbb{R}^n,$ the so-called "column space" of $X.$ Linear algebra tells us the dimension of this space equals the dimension of the domain $(k)$ minus the dimension of the kernel (the "null space" of $X$), whence the dimension cannot exceed $k$ and equals $k$ precisely when the only vector $\beta$ for which $X\beta=0$ is the zero vector itself.

The dimension of the space

Can more be said? Yes, but it's a little delicate. Let's focus on the columns generated by monomials. Suppose they correspond to $x^{e_1}, x^{e_2},$ and so on, through $x^{e_p},$ where all the exponents $e_j$ are distinct positive integers. We may order the exponents so that $1\le e_1 \lt e_2 \lt \cdots \lt e_p,$ whence $e_p \ge p.$

Notice that when a value of $x$ is the same for two observations, these monomials generate identical rows, which therefore are linearly dependent: the extra rows do not change the dimension of the space of the model. Thus, we only need to consider the distinct values of $x$ in the dataset. Let there be $m\le n$ of these and let $X_0$ designate the model matrix for the $m$ distinct observations.

(Geometrically, the equality of some of the original values $x$ means the data reside in some $m$-dimensional subspace $V^m$ of $\mathbb{R}^n.$ Limiting to a set of unique rows of $X$ induces a linear isomorphism $\iota: V^m\to\mathbb{R}^m.$ $X_0$ thereby represents a linear transformation $X_0:\mathbb{R}^p \to \mathbb{R}^m$ given by applying $X$ and following that by $\iota.$ You did ask "precisely what vector space" was involved!)

To compute the dimension of the model space, let's do it indirectly by finding some null vectors. That is, we seek $p$-vectors $\beta$ for which $X_0\beta=0.$ Writing these distinct values of $x$ as $x_1,x_2, \ldots, x_m,$ we have

$$X_0 = \left(x_i^{e_j}\right) = \pmatrix{x_1^{e_1} & x_1^{e_2} & \cdots & x_1^{e_p}\\ x_2^{e_1} & x_2^{e_2} & \cdots & x_2^{e_p}\\ \vdots & \vdots & \ddots & \vdots\\ x_m^{e_1} & x_m^{e_2} & \cdots & x_m^{e_p}}.$$

The equation $X_0\beta = 0$ represents the system of simultaneous equations

$$\beta_1 x_i^{e_1} + \beta_2 x_i^{e^2} + \cdots + \beta_p x_i^{e^p} = 0,$$

$i = 1, 2, \ldots, m.$ In other words,

All the $x_i$ are roots of the polynomial $\beta_1 T^{e_1} + \beta_2 T^{e_2} + \cdots + \beta_p T^{e_p}.$

Because the $x_i$ are distinct, this polynomial has at least $m$ roots. The Fundamental Theorem of Algebra, which asserts any polynomial of degree $p$ (with real or complex coefficients) has at most $p$ distinct roots (real or complex), implies $e_p \ge m.$ Consequently,

The kernel of the model matrix is $\{0\}$ whenever $e_p\lt m:$ that is, whenever the maximum power of your monomials is less than the number of distinct values of $x,$ the dimension of the space generated by those monomials equals the number of monomials. There is no redundancy.

An important corollary is that when the powers are $e_j = 1,2,\ldots, p,$ then so long as you use fewer monomials $(p)$ than distinct values of $x$ $(m),$ the dimension of the model space will equal the number of monomials.

As a simple common example, suppose $x$ is a binary indicator, thereby limited to two values, and that both values appear in the data. This result shows that when you are including an intercept in the model (essentially $x^0=1$) along with $x$ itself, it's pointless to introduce any more monomials.

Constructing a model matrix

It's usually not a good idea to work with the matrices $X$ and $X_0$ whose columns are monomials, because they often create numerical instabilities. Part of the problem is that high powers grow very quickly. If you have, say, distinct values of $x$ from $1$ through $10$ and you choose to include a monomial of degree $e_j=16$ or higher, $1^{e_j}$ is indistinguishable from $0$ in comparison to $10^{e_j}$ when you are using double-precision floating point arithmetic, as illustrated in Python 3.0 here:

Input and output of floating point calculations

Standardizing the columns helps a little, but not much, because the monomials tend to be nearly collinear anyway. See Why are there large coefficents for higher-order polynomial for some discussion.

The standard solution is to compute orthogonal polynomials. These are linear combinations of the monomials $T^{e_j}$ that, when applied to the values $x_i,$ give vectors that are orthogonal to each other. (They form a different basis for the model space, but the model space is unchanged.) The orthogonality is the opposite of nearly collinear: these new columns don't interfere with each other in the underlying calculations. For much more about this, including some statistical subtleties, please read Raw or orthogonal polynomial regression?.

An often overlooked point is that despite being called "orthogonal polynomials," which sounds kind of universal, the orthogonal polynomials depend on the specific distinct values $x_i.$ Some of the links I provide (above and below) mention good ways to compute the orthogonal polynomials for a dataset.

We have many good discussions and examples of the use of orthogonal polynomials in regression models: please see the results of this site search for some of the better answers.

Finally, if you would like to see explicit examples of model matrices involving orthogonal polynomials, I included some in a post on multivariate orthogonal polynomials.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.